#### Abstract

The transverse momentum spectra of , , and produced in proton-proton , proton-antiproton , proton-lead , gold-gold , and lead-lead () collisions over a wide energy range are analyzed by the (two-component) Erlang distribution, the Hagedorn function (the inverse power-law), and the Tsallis-Levy function. The initial temperature is obtained from the color string percolation model from the fit by the (two-component) Erlang distribution in the framework of a multisource thermal model. The excitation functions of several parameters such as the mean transverse momentum and initial temperature increase from 39 GeV to 13 TeV, which is considered in this work. The mean transverse momentum and initial temperature decrease (increase slightly or do not change significantly) with the increase of rapidity (centrality). Meanwhile, the mean transverse momentum of is larger than that of and , and the initial temperature for emission is higher than that for and emission, which shows a mass-dependent behavior.

#### 1. Introduction

The excitation functions of some physical quantities are significative to help us to understand the nuclear reaction mechanism and the system evolution characteristic. For instance, the higher the mean transverse momentum is, the higher excitation state the emission source stays at. Meanwhile, the higher the initial temperature () [1–5] is, the more violent the collisions are. By the analysis of the excitation functions of and , we can learn more about the process in high energy collisions in which the excitation functions of several parameters such as and can be obtained from the spectra of produced particles.

In a data-driven reanalysis, to obtain and , at the first place, we need the spectra of particles in experiments. At the second place, we should choose appropriate functions such as the Erlang distribution [6–8], the Hagedorn function or the inverse power-law [9, 10], and the Tsallis-Levy function [11, 12]. At the last place, we use the chosen functions to fit the experiential data on particle spectra. By describing the spectra, the parameters from the selected functions can be extracted. By comparing the parameters obtained from the experiential data at different energies, centralities, and rapidities, we can find out the dependences of parameters on these quantities. These dependences are related to excitation and expansion degrees of emission source, which is beneficial for us to understand the mechanism and characteristic of nuclear reactions and system evolution.

Besides the two derived parameters and , we can obtain other related parameters by using the method which is similar to extract and . For example, using the Hagedorn function or the inverse power-law [9, 10] and the Tsallis-Levy function [11, 12] to fit spectra, some free parameters such as , , , and in the mentioned functions which will be discussed in Section 2 can be extracted. These free parameters are also useful to understand particle productions and system evolution. Not only the excitation functions of derived parameters and but also the trends of free parameters , , , and can be studied from the fit to spectra.

In this work, the (two-component) Erlang distribution [6–8], Hagedorn function (the inverse power-law) [9, 10], and Tsallis-Levy function [11, 12] are introduced firstly in Section 2. Then, in Section 3, the three distributions or functions are used to preliminarily fit the spectra of heavy flavor quarkonia (charmonia and bottomonia) produced in high energy collisions. The function results are compared with the spectra of , , and measured by the STAR [13, 14], CDF [15–17], ALICE [18], LHCb [19–27], ATLAS [28–30], and CMS Collaborations [31–34] over a wide energy range. Finally, in Section 4, we give our summary and conclusions.

#### 2. Formalism and Method

##### 2.1. The (Two-Component) Erlang Distribution

According to the multisource thermal model [6–8], a given particle is produced in the collision process where a few partons or quarks have taken part in. Each (the -th) parton is assumed to contribute to an exponential function of transverse momentum () distribution. Let denotes the mean transverse momentum contributed by the -th parton, we have the probability density function of to be which is normalized to 1. The probability density function of contributed by all partons which have taken part in the collision process is the convolution of exponential functions [6–8]. We have the distribution (the probability density function of ) of final state particles to be the Erlang distribution

which is naturally normalized to 1. The mean is .

In the two-component Erlang distribution, we have

where denotes the contribution fraction of the first component, () denotes the number of partons in the first (second) component, and denotes the mean transverse momentum contributed by each parton in the first (second) component. The mean is , where –3 in this work and if .

##### 2.2. The (Two-Component) Hagedorn Function

The Hagedorn function is an inverse power-law which is suitable to describe wide spectra of particles produced in the hard scattering process. In refs. [9, 10], the Hagedorn function or the inverse power-law shows the probability density function of to be

where and are the free parameters and is the normalization constant which is related to and and results in . Equation (6) is an empirical formula inspired by quantum chromodynamics (QCD). We call Eq. (4) the Hagedorn function or the inverse power-law [9, 10].

In the case of using two-component Hagedorn function, we have where denotes the contribution fraction of the first component, () is the normalization constant which results in the first (second) component to be normalized to 1, and () and () are free parameters related to the first (second) component. To combine the free parameters of the two components, we have and .

Generally, Eq. (4) is possible to describe the spectra in both the low- and high- regions. In fact, the spectra in the low- and high- regions represent similar trend in some cases. This is caused due to the similarity [35–45] which is widely existent in high energy collisions, where the similarity means the common or universality laws existed in different processes or collisions. In addition, one can revise Eq. (4) if needed in different ways [46–52] which suppress in the spectrum itself in low- or high- region according to the experimental spectra. To discuss various revisions of the Hagedorn function or the inverse power-law [9, 10] is beyond the focus of this paper. We shall not discuss anymore on this issue. For a very wide spectrum, Eq. (5) is possibly needed.

##### 2.3. The (Two-Component) Tsallis-Levy Function

The Tsallis statistics [11] has wide applications in high energy collisions. There are various forms of the Tsallis distribution or function. In this work, we use the Tsallis-Levy function [12]. where and are free parameters, is the transverse mass, is the rest mass of the considered particle, and is the normalized constant which is related to , , and and results in .

We notice that is related to particle mass , which is not the case of and presented in Eqs. (2) and (4), respectively. Although is related to , this relation is not strong due to appearing only in . The fact that the Tsallis distribution depends on shows that this takes simple kinematics into account, as it is well known that or (something like transverse kinetic energy) is a better “scaling variable” for the spectra than .

In the case of using two-component Tsallis-Levy function, we have where denotes the contribution fraction of the first component, () is the normalization constant which results in the first (second) component to be normalized to 1, and () and () are free parameters. To combine the free parameters of the two components, we have and .

The temperature parameter in the Tsallis-Levy function is an effective temperature at the final state (the stage of kinetic freeze-out). This effective temperature is not a “real” temperature because it includes not only the contribution of random thermal motion but also the contribution of flow effect. In the case of the first (second) component having () with the fraction of (), the common effective temperature of the two components is extracted from the assumed common equilibrium state of the two components. That is which has the same form as the parameter .

##### 2.4. The Initial Temperature

According to the color string percolation model [53–55], the initial temperature of the emission source is determined by where is the square of the root-mean-square of due to . If the -component () and -component () of the transverse momentum are considered, we have

In the source rest-frame and under the assumption of isotropic emission, if the -component of momentum is , we also have

Although the source rest-frame is the lab-frame for symmetric collisions, we have mentioned the source rest-frame because asymmetric proton-lead (+Pb) collisions are also considered in this work.

It should be noted that we have used a single string in the cluster for a given particle production because only a projectile participant quark and a target participant quark are mainly considered in our treatment. The assumption of the single string results in the color suppression factor to be 1 in the color string percolation model [54]. If we consider more than one strings taking part in the given particle production, the minimum will be nearly 0.6 [54]. Thus, we shall obtain a higher by multiplying a revised factor in Eqs. (8), (10), and (11). In our opinion, although more than one strings have influences on the given particle production, the main role is the single string.

##### 2.5. Discussion on the Functions

We would like to point out that the three types of functions are mainly just used here as parametrizations to achieve a good fit to the data, to be able to extract and , though the Hagedorn and Tsallis-Levy functions are physically relevant. In fact, in the two functions, if we let , , , the two functions are the same. Here, is an entropy index that characterizes the excitation degree of the collision system [11, 12]. Generally, or is a sizeable quantity, which results in to close to 1 and the collision system to close to an equilibrium state.

We have used the two-component functions in some cases. The reason for using two-component source, i.e., basically two temperatures is not just used to achieve a better fit to the data. Physically, the first component corresponds to the non-head-on collisions between projectile and target participant quarks. The second component corresponds to the head-on collisions between the two quarks. Generally, the first component has a large fraction and low and . The second component has a less fraction and high and . Because the head-on collisions between the two quarks are infrequent, single component function is usually applicable.

In principal, no matter what functions are used to fit the experimental data, (or ) obtained from different fits is approximately the same within a small systematic uncertainty, if different functions fit the data good enough in the region of data available. For example, if simple Maxwell-Boltzmann or Bose-Einstein statistics can fit the data, we may obtain similar (or ) with other functions. In the case of multicomponent Maxwell-Boltzmann or Bose-Einstein statistics being needed, we may also obtain similar (or ).

Indeed, the data itself decides (or ), and (or ) can be directly obtained from the data itself. The reason why we use functions is to see the tendency where the data is not available. However, the extrapolation on the tendency should be careful because it is not fully true, as it to the low- and high regions (where there is no data) could in principle have a major effect on the tendency (for example in case of very step exponentials near , or power-law tails at large ). To reduce the effect, the data should be measured in a sufficiently large interval so that the extrapolation does not spoil as far as possible.

For different components and functions, we do not need to consider the values of mid-rapidity (mid-) or mid-pseudorapidity (mid-), or the values of mid- or mid- can be regarded as 0 directly. In fact, for the experimental data with non-zero mid- or mid-, we may directly regard them as those with mid- or mid-. This treatment is performed to subtract the contribution of kinetic energy of directed motion to the temperature.

#### 3. Results and Discussion

Ordered by center-of-mass energy per nucleon pair or if only one pair) for different panels, Figure 1 shows the spectra, (a–c) , (d, e) , (f) , and (g) , of (a–d) , (e) , (f) prompt , and (g) inclusive produced in (a–c) gold-gold (), (d, e) proton-proton (+), (f) proton-antiproton (), and (g) lead-lead () collisions at mid-rapidity (a–d) , (e) , (f) , and forward rapidity (g) at or (a) 39, (b) 62.4, (c) 200, (d) 500, and (e) 510 GeV, as well as (f) 1.8 and 1.96, and (g) 2.76 TeV, where denotes the number of particles, denotes the cross section, and denotes the yield. The symbols represent the experimental data [13–16, 18] and the curves are our fitted results. In the calculations, the method of least square is used to obtain the best free parameters. The values of free parameters , , , and are listed in Table 1 with and number of degrees of freedom (ndof). The values of free parameters , , , and are listed in Table 2 with and ndof. One can see that the (two-component) Erlang distribution, the Hagedorn function, and the Tsallis-Levy function fit approximately the experimental spectra of via different decay or production modes in high energy , , , and collisions.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

The spectra, , of (a, c, e, and g) prompt and (b, d, f, and h) from produced in (a, b) and (c–h) collisions at or (a, b) 5, (c, d) 7, (e, f) 8, and (g, h) 13 TeV are presented in Figure 2. The symbols represent the experimental data [19–22], and the curves are our fitted results. The method of least square is used to obtain the best parameter values which are listed in Tables 1 and 2 with and ndof. One can see that the experimental spectra of via different production modes in different rapidity intervals in and collisions at high energies are approximately fitted by the Erlang distribution, the Hagedorn function, and the Tsallis-Levy function.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

Figure 3 shows the spectra, (a, c–g) and (b) , of via different production modes. The symbols represent the experimental data [15, 16, 23, 28, 31, 32], and the curves are our fitted results. The values of free parameters are listed in Tables 1 and 2 with and ndof. One can see that the experimental spectra of via different production modes in different rapidity intervals in , , and collisions at high energies are approximately fitted by the (two-component) Erlang distribution, the Hagedorn function, and the Tsallis-Levy function.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

Figure 4 shows the spectra, , of via different production modes in collisions at (a, b) 7, (c) 8, and (d) 13 TeV. The symbols represent the experimental data [29, 30, 33], and the curves are our fitted results. The values of free parameters are listed in Tables 1 and 2 with and ndof. For Figure 4(d), the two-component Eqs. (5) and (7) are used, where the free parameters for the first (second) component are listed in the first (second) row. One can see that the experimental spectra of via different production modes in different rapidity intervals in collisions at high energies are also approximately fitted by the Erlang distribution, the (two-component) Hagedorn function, and the (two-component) Tsallis-Levy function.

**(a)**

**(b)**

**(c)**

**(d)**

In Figure 5, the spectra, (a, c–f) and (b) , of (a–c) , (d) , (e) , and (f) induced in (a) and (b–f) collisions at (a) 1.8, (b) 2.76, (c) 5.02, and (d–f) 7 TeV are given. The symbols represent the experimental data [17, 24, 25, 34], and the curves are our fitted results. The parameter values are listed in Tables 1 and 2 with and ndof. One can see that the experimental spectra of in different rapidity intervals in and collisions at high energies are approximately fitted by the Erlang distribution, the Hagedorn function, and the Tsallis-Levy function.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

In Figure 6, the spectra, , of (a, d, and e) , (b, f) , and (c, g) induced in (a–c) and (e–g) + and (d) collisions at or (a–c) 8, (d) 8.16, and (e–g) 13 TeV are given. The symbols represent the experimental data [25–27], and the curves are our fitted results. The parameter values are listed in Tables 1 and 2 with and ndof. Once again, one can see that the experimental spectra of in different rapidity intervals in and collisions at high energies are approximately fitted by the Erlang distribution, the Hagedorn function, and the Tsallis-Levy function.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

Before discussing the trends of parameters, we would like to point out the usability of the concept of temperature in () collisions, which are small in size. As in refs. [56–59], in this work, we have treated + (+) collisions as where a medium was formed, or at least there is some degree of thermalization, enough to have a temperature for the emission source. On the other hand, the temperature parameter of the emission source is a reflection of the average kinetic energy of given particles. This means that we may use the concept of temperature. Even if the collision system is not enough large, we may use the temperature parameter to characterize the average kinetic energy of given particles over many events.

Figure 7 shows the dependences of (a, c, and e) and (b, d, and f) on (or ) for (a, b) , (c, d) , and (e, f) . The different symbols represent the parameter values derived from free parameters extracted from Figures 1–6 and listed in Tables 1 and 2, where only the (two-component) Erlang distribution in the region of data available is used as an example. It is expected that the results corresponding to the Hagedorn and Tsallis-Levy functions are very close to the plot, because the two functions also describe approximately the data. As what we discussed in Subsection 2.5, no matter what functions are used to fit the experimental data, one should obtain similar (or ), if different functions fit the data good enough in the region of data available. By using the mentioned three functions which fit the data good enough, one can obtain (or ) within a systematic uncertainty of 8%. One can see from Figure 7 that and increase significantly with the increase of collision energy. Meanwhile, and increase with the increase of particle mass.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Figure 8 is the same as Figure 7, but showing the dependences of (a, c, e, and g) and (b, d, f, and h) on (a, b) centrality and (c-h) rapidity for (a–d) , (e, f) , and (g, h) . The different symbols represent the parameter values derived from free parameters extracted from Figures 1–6 and listed in Tables 1 and 2, where only the (two-component) Erlang distribution is used as an example. One can see that and increase slightly with the increase of event centrality from peripheral to central collisions and decrease with the increase of rapidity from mid-rapidity to forward rapidity. Meanwhile, and increase with the increases of collision energy and particle mass.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

The above parameter tendencies show that the temperature is mass dependent. This is also a reflection of the formation time dependence. According to the hydrodynamic behavior, “massive particles coming out of the system earlier in time with smaller radial flow velocities” [59]. This means that with the increase of mass, the formation time decreases, the temperature increases, and the flow velocity decreases. It should be noted that the fact that the massive particles coming out the system earlier is not caused by the high excitation of the system, but the leaver over due to the inertia of massive particles, in the hydrodynamic evolution.

We may explain the tendency of derived and which have similar tendency with . With the increase of collision energy, the violent degree of collisions increases significantly due to large energy transfer, which results in the obvious increase of . With the increase of centrality, the degree of multiple-scattering increases due to more participant nucleons and produced particles taking part in the scattering process, which results in a slight increase of emission angle and then a slight increase of . With the increase of rapidity, the energy transfer decreases due to larger penetrability between participant nucleons. Meanwhile, the degree of multiple-scattering also decreases due to less produced particles taking part in the scattering process. These two factors result in the decrease of . It is natural that increases with the increase of .

The free parameters in the Erlang distribution are directly reflected in , which will not be discussed anymore. The free parameters in the Hagedorn and Tsallis-Levy functions will be discussed in the Appendix, because no clear physical conclusion can be drawn from them at present, though the tendencies of them can be seen from this work.

#### 4. Summary and Conclusions

In summary, the transverse momentum spectra of , , and produced in , , , , and collisions over an energy range from 39 GeV to 13 TeV have been analyzed by the (two-component) Erlang distribution, the Hagedorn function, and the Tsallis-Levy function. The function results are approximately in agreement with the experimental data measured by several international collaborations. The values of related parameters are extracted from the fits, and the excitation functions of these parameters are obtained.

The excitation functions of parameters and increase from 39 GeV to 13 TeV. Meanwhile, and increase (slightly) with event centrality and particle mass and decrease from mid-rapidity to forward rapidity. These tendencies render that these parameters describe the excitation and expansion degrees of the system. At higher energy, larger energy transfer had happened, which results in higher excitation and expansion degrees of the system. In central collisions and at mid-rapidity, larger energy transfer and further multiple-scattering had happened, which also results in higher excitation and expansion degrees of the system.

The parameters () and () increase with the collision energy, which reflects the degree of energy deposition and transfer. In given collisions, there is a negative correlation between () and (). At different energies, there is a positive correlation between () and (). Indeed, there are correlations between () and () when we determine these parameters. The correlation between () and () is similar to that between kinetic freeze-out temperature and transverse flow velocity. If () is similar to kinetic freeze-out temperature, () should be similar to transverse flow velocity.

#### Appendix

#### A. Figures of Parameters and Some Discussions

Figure 9 is the same as Figure 7, but showing the dependences of (a, c, e) and (b, d, and f) on (or ) for (a, b) , (c, d) , and (e, f) . The different symbols represent the parameter values derived from free parameters extracted from Figures 1–6 and listed in Tables 1 and 2. One can see that and increase with the increase of collision energy and particle mass.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Figure 10 is the same as Figure 7, but showing the dependences of (a, c, e, and g) and (b, d, f, and h) on (a, b) centrality and (c–h) rapidity for (a–d) , (e, f) , and (g, h) . One can see that increases (decreases) slightly with the increase of event centrality from peripheral to central collisions and decreases (increases) (slightly) with the increase of rapidity from mid-rapidity to forward rapidity. Meanwhile, increases with the increases of collision energy and particle mass.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

Figure 11 is the same as Figure 7, but showing the dependences of (a, c, and e) and (b, d, and f) on (or ) for (a, b) , (c, d) , and (e, f) . One can see that and increases with the increase of collision energy and particle mass.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Figure 12 is the same as Figure 7, but showing the dependences of (a, c, e, and g) and (b, d, f, and h) on (a, b) centrality and (c–h) rapidity for (a–d) , (e, f) , and (g, h) . One can see that increases (decreases) slightly with the increase of event centrality from peripheral to central collisions and decreases (increases) (slightly) with the increase of rapidity from mid-rapidity to forward rapidity. Meanwhile, increases with the increases of collision energy and particle mass.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

**(g)**

**(h)**

The tendency of and with collision energy are also explained by more violent collision at higher energy. Both and increase with the increase of collision energy. This means that and are positively correlative at different energies. Meanwhile, for a given spectrum or in given collisions, an increase in is concomitant with a decrease in . This means that and are negatively correlative in given collisions (at given energy). There are correlations between and when we determine these parameters.

The correlation between and is similar to that between kinetic freeze-out temperature and transverse flow velocity [60, 61] which also show positive correlation at different energies and negative correlation in a given spectrum. If is similar to kinetic freeze-out temperature, should be similar to transverse flow velocity. Meanwhile, the results obtained in this work are in agreement with our recent work [62], which shows mass-dependent parameters. In particular, with the increase of particle mass, , , , and increase.

#### Data Availability

The data used to support the findings of this study are included within the article and are cited at relevant places within the text as references.

#### Ethical Approval

The authors declare that they are in compliance with ethical standards regarding the content of this paper.

#### Disclosure

The funding agencies have no role in the design of the study; in the collection, analysis, or interpretation of the data; in the writing of the manuscript; or in the decision to publish the results.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant Nos. 11575103 and 11947418, the Scientific and Technological Innovation Programs of Higher Education Institutions in Shanxi (STIP) under Grant No. 201802017, the Shanxi Provincial Natural Science Foundation under Grant No. 201901D111043, and the Fund for Shanxi “1331 Project” Key Subjects Construction.