Excitation Functions of Tsallis-Like Parameters in High-Energy Nucleus–Nucleus Collisions

The transverse momentum spectra of charged pions, kaons, and protons produced at mid-rapidity in central nucleus–nucleus (AA) collisions at high energies are analyzed by considering particles to be created from two participant partons, which are assumed to be contributors from the collision system. Each participant (contributor) parton is assumed to contribute to the transverse momentum by a Tsallis-like function. The contributions of the two participant partons are regarded as the two components of transverse momentum of the identified particle. The experimental data measured in high-energy AA collisions by international collaborations are studied. The excitation functions of kinetic freeze-out temperature and transverse flow velocity are extracted. The two parameters increase quickly from ≈3 to ≈10 GeV (exactly from 2.7 to 7.7 GeV) and then slowly at above 10 GeV with the increase of collision energy. In particular, there is a plateau from near 10 GeV to 200 GeV in the excitation function of kinetic freeze-out temperature.


Introduction
High-energy collider experiments are designed to study the strongly interacting matter at high temperatures and densities [1]. The deconfinement of colliding hadrons into quarkgluon plasma (QGP), which then rapidly expands and cools down [2], is conjectured to be created at such extreme collision energies [3][4][5][6]. In high-energy and nuclear physics, the study of transverse (momentum (p T ) or mass (m T )) spectra of charged particles produced in nucleus-nucleus (AA) collisions is very important. In particular, the AA collision process at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron Collider (LHC) provides a good opportunity to study the signals and characteristics of QGP generation, so as to indirectly study the system evolution and the reaction mechanism of particle generation.
During the time evolution of collision system [7][8][9], the stages of kinetic freeze-out and chemical freeze-out are two important processes. In the stage of chemical freeze-out, a phase transition from QGP to hadrons occurred in the system, so the composition and ratio of various particles remain unchanged. In the stage of kinetic freeze-out, elastic collisions among particles stop, so their p T and then m T spectra are unchanged [8,10]. Therefore, by studying the p T (m T ) spectra, we can obtain some useful information, such as the effective temperature (T), the chemical freeze-out temperature (T ch ), and the kinetic freeze-out temperature (T 0 or T kin ) of the system, as well as the transverse flow velocity (β T ) of the final state particles. The temperature in which we do not exclude the contribution of transverse flow is called the effective temperature, which is related to the kinetic freeze-out the experimental data measured by the E866 [31], E895 [32,33], E802 [34,35], NA49 [36,37], STAR [38][39][40], and ALICE Collaborations [41][42][43], we analyze the tendency of parameters.
The remainder of this paper is structured as follows. The formalism and method are shortly described in Section 2. Results and discussion are given in Section 3. In Section 4, we summarize our main observations and conclusions.

Formalism and Method
The Tsallis distribution has different forms or revisions [44][45][46][47], we have the Tsallis-like distribution of p T at mid-y to be where N denotes the number of particles, can be obtained using p T , is an entropy index that characterizes the degree of equilibrium or non-equilibrium, n is a parameter related to q, and µ is the chemical potential. In particular, in the expression of m T − µ − m 0 , m T is simplified from m T cosh y because cosh y ≈ 1 at mid-y.
We have the probability density function of p T at mid-y to be 1 N dN dp T ∝ m T . (4) Empirically, to fit the spectra of p T at mid-y in this work, Equation (4) can be revised as where C is the normalization constant; a 0 is a new non-dimensional parameter that describes the bending degree of the distribution in low-p T region (p T = 0 ∼ 1 GeV/c), which is introduced artificially and tested in our recent work [48,49]; and m a 0 T is revised from m T due to the introduction of the revised index a 0 . Because of the limitation of the normalization, changing the bending degree in the low-p T region will change the slope in the high-p T region. Although writing Cm a 0 T in Equation (5) is not ideal, as it yields a fractional power unit in C, we have no suitable method to scale out the unit by e.g., m 0 due to the nonlinear relationship between m T and m 0 shown in Equation (2). In Equation (5), the other parameters such as q and a 0 do not appear in the function name for the purpose of convenience. In this work, we call Equation (5) the revised Tsallis-like function.
In the framework of the multisource thermal model [30], we assume that two participant partons take part in the collisions. Let p t1 and p t2 denote the components contributed by the first and second participant parton to p T , respectively, where p t1 (p t2 ) is less than the transverse momentum of the participant parton. We have where the two components are perpendicular due to the fact that p t1 and p t2 are assumed to be the two components of the vector p T . Although multiparton collisions can be important especially for central high-energy nucleus-nucleus collisions, the main contributors to particle production are still binary parton collisions, which are also the basic collision process. After all, the probability that three or more partons collide simultaneously is small. Instead, the probability of binary parton collisions is large. In binary parton collisions, each parton, e.g., the i-th parton, is assumed to contribute to p T to obey Equation (5), where i = 1 and 2. The probability density functions at mid-y obeyed by p t1 and p t2 is where the subscript i is used for the quantities related to the i-th parton and m 0i is empirically the constituent mass of the considered parton. Generally, in the case of considering u and/or d quarks, we take m u = m d = 0.3 GeV/c 2 . It is noted that the constituent quark masses of 0.3 GeV are not incompatible with the pion and kaon masses because the collisions between the two participant quarks can produce more than one particle. The conservation of energy is satisfied in the collisions. The value of µ i will be discussed at the end of this section. Let φ denote the azimuthal angle of p T relative to p t1 . According to the works in [50,51], we have the unit normalized probability density function of p T and φ to be where f 1,2 (p t1 , p t2 , T) denotes the united probability density function of p t1 and p t2 . Further, we have the probability density function of p T to be Equation (9) can be used to fit the p T spectra and obtain the parameters T, q, and a 0 . In the case of fitting a wide p T spectrum e.g., p T > 5 GeV/c, Equation (9) cannot fit well the spectra in high-p T region. Then, we need a superposition of one Equation (9) with low T and another Equation (9) with high T to fit the whole p T spectrum. As will be seen in Figure 3e in the next section, the contribution fraction of the low T component is very large (≈99.9%). In most cases in Figures 1-3, we do not need the superposition due to narrow p T spectra. In the case of using a two-component distribution, we have the probability density function of p T to be where k (1 − k) denotes the contribution fraction of the first (second) component and ] is given by Equation (9). The second component is related to the core-corona picture as mentioned later on in detail in subsection 3.3. Correspondingly, the temperature is averaged by weighting the two fractions. The temperature T defined by Equation (11) reflects the common effective temperature of the two components which are assumed to stay in a new equilibrium in which T still characterizes the average kinetic energy. Similarly, the weighted average can be performed for other parameters in the two components in Equation (10). Note that the limit of the first and second (low-and high-p T ) components is determined by a convenient treatment. Generally, the contribution fraction k of the first component should be taken as largely as possible. As will be seen in the next section, we take k = 1 in most cases; only in Figure 3e we take k = 0.999. Because the contribution fraction of the second component is zero or small enough, Equation (10) becomes Equation (9), and the weighted average of the two parameters in Equation (10) becomes the parameter in Equation (9). Because Equations (1), (4), (5), and (7) are suitable at mid-y, Equations (8)- (10) are also suitable at mid-y. In addition, the rapidity ranges quoted in the next section are narrow and around 0, though the concrete ranges are different. This means that the mentioned equations are applicable.
We would like to point out that although the model used by itself is not enough to provide information of the deconfinement phase transition from hadronic matter to QGP, the excitation function of extracted parameter is expected to show some particular tendencies. These particular tendencies include, but are not limited to, the peak and valley structures, the fast and slow variations, the positive and negative changes, etc. These particular tendencies are related to the equation of state (EOS) of the considered matter. The change of EOS reflects the possible change of interaction mechanism from hadron-dominated to parton-dominated intermediate state. Then, the deconfinement phase transition of the considered matter from hadronic matter to QGP is possibly related to the particular tendencies. It is natural that the explanations are not only for a given set of data. The present model will show a method to fit and explain the data.
To obtain β T , we need to know the slope of p T versus m in the source rest frame of the considered particle. That is, we need to calculate p T and m. According to Equation (10), we have due to where p T max denotes the maximum p T .
As the mean energy, E = m = p 2 + m 2 0 , where p is the momentum of the considered particle in the source rest frame. The analytical calculation of m is complex. Instead, we can perform the calculation by the Monte Carlo method. Let R 1,2 denote random numbers distributed evenly in [0, 1]. Each concrete p T satisfies where δp T denotes a small shift relative to p T . Each concrete emission angle θ satisfies θ = 2 arcsin R 2 (15) due to the fact that the particle is assumed to be emitted isotropically in the source rest frame. Each concrete momentum p and energy E can be obtained by and respectively. After repeating calculations multiple times in the Monte Carlo method, we can obtain E, that is, m. Then, the slope of p T versus m is identified as β T . Meanwhile, the intercept of T versus m 0 is identified as T 0 . Here, we emphasize that we have used the alternative method introduced in Section 1 to obtain T 0 and β T .
Note that in some cases, the transverse spectra are shown in terms of m T , but not p T . To transform the probability density function f p T (p T , T) of p T to the probability density function f m T (m T , T) of m T , we have the relation Then, we have due to Equation (2). Using the parameters from m T spectra, we may also obtain T 0 , p T , m, and β T . We now discuss the chemical potential µ i of the i-th parton. Generally, the chemical potential µ of a particle obviously affects the particle production at low energy [52][53][54][55][56][57][58].
For baryons (mostly protons and neutrons), the chemical potential µ B related to collision energy √ s NN is empirically given by where both µ B and √ s NN are in the units of GeV [59][60][61]. According to the authors of [52], we have µ u = µ d = µ B /3 because a proton or neutron is consists of three u/d quarks (i.e., uud or udd).

Comparison with Data and Tendencies of Free Parameters
, of charged pions, kaons, and protons produced in 0-5% Au-Au, Pb-Pb, and Xe-Xe collisions at different √ s NN . The collision types, particle types, mid-y ranges, centrality classes, and √ s NN are marked in the panels. The symbols represent the experimental data measured by different collaborations. The solid and dashed curves are our results, fitted by using Equation (10) due to Equations (7) and (9), with µ i = 0 and µ i = µ B /3, respectively. In the process of fitting the data, we determine the best parameters by the method of least squares. The experimental uncertainties used in calculating the χ 2 are obtained by the root sum square of the statistical uncertainties and the systematic uncertainties. The parameters that minimize the χ 2 are the best parameters. The errors of parameters are obtained by the statistical simulation method [62,63], which uses the same algorithm as usual, if not the same Code, in which the errors are also extracted from variations of χ 2 . The values of T 1 , T 2 , k, q, and a 0 are listed in Tables 1 and 2 with the normalization constant (N 0 ), χ 2 , and the number of degree of freedom (ndof), or explained in the caption of Table 1.
In a few cases, the values of χ 2 /ndof are very large (5-10 or above), which means "bad" fit to the data. In most cases, the fits are good due to small χ 2 /ndof which is around 1. To avoid possible wrong interpretation with this result, the number of "bad" fits are limited to be much smaller than that of good fits, for example, 1 to 5 or more strict such as 1 to 10. Meanwhile, we should also use other method to check the quality of fits. In fact, we have also calculated the p-values in the Pearson method. It is shown that all p-values are less than 3 × 10 −7 . These p-values corresponds approximately to the Bayes factor being above 100 and to the confidence degree of 99.99994% at around 5 standard deviations (5σ) of the statistical significance. This means that the model function is in agreement with the data very well. To say the least, most fits are acceptable.
Note that we will use a set of pion, kaon, and proton spectra to extract T 0 and β T in Section 3.2. For energies in the few GeV range, the spectra of some negative particles are not available in the literature. Therefore, we have to give up to analyze all the negative particle spectra in Figure 1. In our recent work [28], the positive and partial negative particle spectra were analyzed by the standard distribution. The tendencies of parameters are approximately independent of isospin, if not the same for different isospins.
In panel (f), the factor 1/2π does not appear, which causes different normalization from other panels. The symbols represent the experimental data at mid-y measured by the E866, E895, and E802 Collaboration at the AGS [31][32][33][34][35] and by the NA49 Collaboration at the SPS [36,37]. The solid and dashed curves are our results, fitted by using Equation (10) due to Equations (7) and (9), with µ i = 0 and µ i = µ B /3, respectively.   included on the vertical axis, which can be omitted. The symbols represent the experimental data at mid-y measured by the STAR Collaboration at the RHIC [38][39][40]. The solid and dashed curves are our results, fitted by using Equation (10) due to Equations (7) and (9), with µ i = 0 and µ i = µ B /3, respectively.
|y|<0.1   In panels (c,d,f), the factor 1/N EV is included on the vertical axis, which can be omitted. In panels (e,f), the item (2π p T ) −1 is not included on the vertical axis, which results in different calculation for vertical values from other panels in the normalization. The symbols represent the experimental data at mid-y measured by the STAR Collaboration at the RHIC [38][39][40] and by the ALICE Collaboration at the LHC [41][42][43]. The solid and dashed curves are our results, fitted by using Equation (10) due to Equations (7) and (9), with µ i = 0 and µ i = µ B /3, respectively. (10) approximately describes the considered experimental data. For all energies and particles, T 1 and T 2 are identical except for the 5.02 TeV Pb-Pb data from ALICE. This means that none of the spectra have a wide enough range to determine the second component except the data at 5.02 TeV. The two-component fit is only really used at 5.02 TeV. In the high-p T region, the hard scattering process which is described by the second component in Equation (10) contributes totally. However, in the case of using the two-component function, k (= 0.999) is very close to 1, which implies that the contribution of the second component is negligible.

One can see from Figures 1-3 and Tables 1 and 2 that Equation
In fact, the second component contributes to the spectrum in high-p T region with small fraction, which does not affect significantly the extraction of parameters. Instead, the parameters are determined mainly by the spectrum in low-p T region. Table 1. Values of free parameters (T 1 , T 2 , q, and a 0 ), normalization constant (N 0 ), χ 2 , and ndof corresponding to the solid curves in Figures 1-3 in which the data are measured in special conditions (mid-y ranges and energies) by different collaborations, where T 2 is not available in most cases because k = 1. In a few cases (at √ s NN = 5.02 TeV), T 2 is available in the next line, where k = 0.999 ± 0.001, which is not listed in the table.  Although the contribution fraction of the second component is very small, the spectra with a wide p T range on Figure 3e are well fit using the two components, this means increasing the number of parameters compared with just Tsallis function. Generally, the spectrum shapes of different particles are different. However, we may use the same function with different parameters and normalization constants to fit them uniformly. In some cases, the spectrum forms are different. We need to consider corresponding normalization treatments so that the fitting function and the data are compatible and concordant.
The value of µ i affects mainly the parameters at below dozens of GeV. Although µ i = 0 is not justified at lower energies, we present the results with µ i = 0 for comparison with µ i = µ B /3 so that we can have a quantitative understanding on the influence of µ i . Note that µ i is only for µ u and µ d , that is µ u = µ d = µ B /3. For pions, we have µ π = µ u + µ d = 2µ B /3. For kaons, we have no suitable expression because the chemical potential µ s for s quark is not available here. Generally, µ s > µ u . Therefore, As a function with wide applications, the Tsallis distribution can describe in fact the spectra presented in Figures 1-3 in most cases, though the values of parameters may be changed. However, to extract some information at the parton level, we have regarded the revised Tsallis-like function (Equation (7)) as the components of p T contributed by the participant partons. The value of p T is then taken to be the root sum square of the components. In the present work, we have considered two participant partons and two components. This treatment can be extended to three and more participant partons and their components. In the case of the analytical expression for more components becoming difficult, we may use the Monte Carlo method to obtain the components, and p T is also the root sum square of the components. Then, the distribution of p T is obtained by the statistical method.
To study the changing trends of the free parameters, Figure 4 shows the dependences of (a) effective temperature T, (b) entropy index q, and (c) revised index a 0 on collision energy √ s NN , where the closed and open symbols are cited from Tables 1 and 2 which are obtained from the fittings with µ i = 0 (solid curves) and µ i = µ B /3 (dashed curves) in Figures 1-3, respectively. The triangles, circles, squares, and pentagrams represent the results for charged pions, kaons, protons, and the average by weighting different yields, respectively. Because the errors of parameters are very small, the error bars in the plots are invisible. One can see from Figure 4 that T increases significantly, q increases slowly, and a 0 increase quickly from ≈ 3 to ≈ 10 GeV (exactly from 2.7 to 7.7 GeV) and then changes slowly at above 10 GeV, except for a large increase (≈ 50%) at the maximum energy, with the increase of ln( √ s NN ). These parameters also show their dependences on particle mass m 0 : With the increase of m 0 , T and a 0 increase and q decreases significantly. Indeed, µ i affects only the parameters at the lower energies (below dozens of GeV), but not higher energy. The behavior of excitation function of T will be discussed as that of T 0 in the next subsection. The large fluctuations of q for pions are caused by the large influence of strong decay from high-mass resonance and weak decay from heavy flavor hadrons. For light particles such as pions, the influence and then the fluctuations are large; while for relative heavy particles such as kaons and protons, the influence and then the fluctuations are small. No matter how large the fluctuations are, the values of q are close to 1.
As we mentioned in the above section, the entropy index q reflects the degree of equilibrium or non-equilibrium of collision system. Usually, q = 1 corresponds to an ideal equilibrium state and q 1 means a non-equilibrium state. The present work shows that q is very close to 1 which means that the system stays in the equilibrium state. Generally, the equilibrium is relative. For the case of non-equilibrium, we may use the concept of local equilibrium. If q is not too large, for example, q ≤ 1.25 or n ≥ 4, the collision system is still in equilibrium or local equilibrium [45,64]. In particular, the system is closer to the equilibrium when it emits protons at lower energy, comparing with pions and kaons at higher energy. The reason is that most protons came from the participant nuclei directly. They have enough time to reach to the equilibrium in the evolution. At lower energy, the system is closer to the equilibrium because the evolution is slower and the system has more time to result in the equilibrium. From the initial collisions to kinetic freeze-out, the evolution time is very short. The lower the collision energy is, the longer the evolution time is. The values of a 0 for the spectra of charged pions, kaons, and protons at above 10 GeV are approximately 0.75, 1, and 1.8, respectively, which drop obviously for pions and kaons at lower energy due to the hadronic phase. In addition, due to the existence of participant protons in both the hadronic and QGP phases, the energy dependence of a 0 for protons is not obvious. Although it is hard to explain exactly the physical meaning of a 0 , we emphasize here that it shows the bending degree of the spectrum in low-p T region [48,49] and affects the slopes in high-p T region due to the limitation of normalization. A large bending degree means a large slope change. In fact, a 0 is empirically related to the contributions of strong decay from high-mass resonance and weak decay from heavy flavor hadrons. This is because that a 0 affects mainly the spectra in low-p T region which is just the main contribution region of the two decays.
One can see that the values of q and a 0 change drastically with particle species. This is an evidence of mass-dependent differential kinetic freeze-out scenario [26]. The massive particles emit earlier than light particles in the system evolution. The earlier emission is caused due to the fact that the massive particles are left behind in the evolution process, but not their quicker thermal and flow motion. In fact, the massive particles have no quicker thermal and flow motion due to larger mass. Instead, light particles have quicker thermal and flow motion and longer evolution time. Finally, light particles reach larger space at the stage of kinetic freeze-out.
The influence of µ i on q and a 0 is very small. Although the prefactor a 0 can come from the Cooper-Frye term (and/or a kind of saddlepoint integration) as discussed, e.g., in [65,66], it is a fit parameter in this work. As an average over pions, kaons, and protons, a 0 is nearly independent of √ s NN at above 10 GeV. As √ s NN increasing from ≈3 to ≈10 GeV, the increase of a 0 shows different collision mechanisms comparing with that at above 10 GeV. Our recent work [67] shows that the energy ≈10 GeV discussed above is exactly 7.7 GeV.

Derived Parameters and Their Tendencies
As we know, the effective temperature T contains the contributions of the thermal motions and flow effect [68]. The thermal motion can be described by the kinetic freeze-out temperature T 0 , and the flow effect can be described by the transverse flow velocity β T . To obtain the values of T 0 and β T , we analyze the values of T presented in Tables 1 and 2, and calculate p T and m based on the values of parameters listed in Tables 1 and 2. In the calculation performed from p T to p T and m by the Monte Carlo method, as in [24][25][26], an isotropic assumption in the rest frame of emission source is used. Figure 5a-f shows the relationship of T and m 0 , determined fitting AA collision systems by our model. Figure 6a Tables 3 and 4. One can see that, in most cases, the mentioned relations are described by a linear function. In particular, the intercepts in Figure 5a-f are regarded as T 0 , and the slopes in Figure 6a-f are regarded as β T , as what we discussed above in the alternative method. Because different "thermometers" are used, T 0 extracted from the intercept exceeds (is not in agreement with) the transition temperature which is independently determined by lattice QCD to be around 155 MeV. To compare the two temperatures, we need a transform equation or relation which is not available at present and we will discuss it later. Table 3. Values of intercepts, slopes, and χ 2 for the solid lines in Figures 5 and 6, where ndof = 1, which is not shown in the table. The units of the intercepts in Figures 5 and 6 are GeV and GeV/c, respectively. The units of the slopes in Figures 5 and 6 are c 2 and c, respectively.

Figure
Relation System √ s N N (GeV) Intercept Slope χ 2 Figure 5a T − m 0    Table 4. Values of intercepts, slopes, and χ 2 for the dashed lines in Figures 5 and 6. It is noted that the above argument on T 0 and β T is based usually on exact hydrodynamic calculations, as, e.g., given in [17,65,[69][70][71][72]. However, in these cases, usually T is extracted, and then some T = T 0 + m 0 u t 2 like correspondence is derived (where instead of m 0 , also energy or average energy could stand, depending on the calculation). Here, as we know, u t is related but not equal to β T , as discussed in the mentioned literature. Therefore, we give up to use u t as β T in this work.
We think that T 0 can be also obtained from p T , and β T can be also obtained from T. However, the relations between T 0 and p T , as well as β T and T, are not clear. Generally, the parameters T 0 and β T are model-dependent. In other models, such as the blast-wave model [17][18][19][20][21], T 0 and β T can be obtained conveniently. The two treatments show similar tendencies of parameters on √ s NN and event centrality, if we also consider the flow effect in small system or peripheral AA collisions [73,74] in the blast-wave model. In order to more clearly see the tendencies of T 0 and β T , we show the dependences of T 0 on √ s NN , β T on √ s NN , and T 0 on β T in Figure 7a-c, respectively. One can see that the two parameters increase quickly from ≈ 3 to ≈ 10 GeV and then slowly at above 10 GeV with the increase of √ s NN in general. There is a plateau from near 10 GeV to 200 GeV.
In particular, T 0 increases with β T due to the fact that both of them increase with √ s NN .
These incremental tendencies show that, in the stage of kinetic freeze-out, the degrees of excitation and expansion of the system increase with increasing √ s NN . These results are partly in agreement with the blast-wave model which shows decreasing tendency for T 0 and increasing tendency for β T with increasing √ s NN from the RHIC [40] to LHC [41] because different partial p T ranges in the data are considered for different particles, while this work uses the p T range as wide as the data. The chemical potential shows obvious influence on T 0 at the lower energies (below dozens of GeV). After considering the chemical potential, the plateau in the excitation function of T 0 becomes more obvious. With the increase of √ s NN , the fact that the values of T 0 and β T increase quickly from ≈ 3 to ≈ 10 GeV and then slowly at above 10 GeV implies that there are different collision mechanisms in the two energy ranges. In AA collisions, if the baryon-dominated effect plays a more important role at below 10 GeV [75], the meson-dominated effect should play a more important role at above 10 GeV. In the baryon-dominated case, less energies are deposited in the system, and then the system has low excitation degree and temperature. In the meson-dominated case, the situation is opposite. Indeed, ≈ 10 GeV is a particular energy which should be paid more attention. It seems that the onset energy of deconfinement phase transition from hadronic matter to QGP is possibly 10 GeV or slightly lower (e.g., 7.7 GeV [67]). If we regard the plateau from near 10 to 200 GeV in the excitation functions of T 0 and β T as a reflection of the formation of QGP liquid drop, the quick increase of T 0 and β T at the LHC is a reflection of higher temperature QGP liquid drop due to larger energy deposition. At the LHC, the higher collision energy should create larger energy density and blast wave, and then higher T 0 and β T . Although any temperature needs to be bound by the phase transition on one side and free streaming on the other side, larger energy deposition at the LHC may heat the system to a higher temperature even the phase transition temperatures at the LHC and RHIC are the same. Both the formed QGP and hadronized products are also possible to be heated to higher temperature.
Although we mentioned that the plateau apparent in T 0 versus √ s NN is possibly connected to the onset of deconfinement, the temperature measured in this work is connected only to T 0 which is usually much smaller than the quark-hadron transition temperature. Because the collision process is very complex, the √ s NN dependence of T 0 reflects only partial properties of the phase structure of a quark medium. To make a determined conclusion, we may connect to the dynamics of the hadron gas. This topic is beyond the focus of the present work and will not be discussed further here. We would like to point out that, in the last three paragraphs mentioned above, the discussions on the excitation function of T 0 presented in Figure 7a are also suitable to the excitation function of T presented in Figure 4a, though the effect of flow is not excluded from Figure 4a. Because the quality of fits is not sufficient in a few cases, our main conclusion that the rise of temperature below 10 GeV suggests that a deconfinement of hadronic matter to QGP is weak. The information of phase transition happened at higher temperatures and near the chemical freeze-out may be reflected at the kinetic freeze-out of a hadronic system. The plateau structure appeared in the excitation function T 0 is expected to relate to the phase transition, though this relation is not clear at present. Other works related to this issue are needed to make a strong conclusion. In other words, to conclude the onset of deconfinement just from the quality of some fits is a loose interpretation. More investigations are needed and also comparison with other findings. This issue is beyond the scope of this analysis.

Further Discussion
The model presented in the analysis can be regarded as a "thermometer" to measure temperatures and other parameters at different energies. Then, the related excitation functions can be obtained, and the differences from the transition around critical point and other energies can be seen. Different models can be regarded as different "thermometers". The temperatures measured by different "thermometers" have to be unified so that one can give a comparison. If we regard the phase transition temperature determined by lattice QCD as the standard one, the values of T 0 obtained in this paper should be revised to fit the standard temperature. However, this revision is not available for us at present due to many uncertain factors. In fact, we try to focus on the "plateau" in the energy dependence of T 0 , but not on the T 0 values themselves.
In addition, the model assumes the contributions from two participant partons in the framework of multisource thermal model [30]. In pp collisions, one can see the point of a hard scattering between two partons and look at the high p T particle productions or other observations. However, even in pp collisions there are underlying events, multiple-parton interactions, etc. Further, the data used in this analysis are from central AA collisions, where hundreds and thousands of hadrons are produced. Although many partons take part in the collisions, only a given two-parton process plays main role in the production of a given set of particles. Many two-parton processes exist in the collisions. Using a model inspired by two participant partons is reasonable.
Of course, one may also expect that the production of many particles can result from three or more partons. If necessary, we may extend the picture of two participant partons to that of three or multiple participant partons [30] if we regard p T of identified particle as the root sum square of the transverse momenta of three or multiple participant partons. It is just that the picture of two participant partons is enough for the production of single particle in this analysis. Besides, we did not try to distinguish between local thermalization of a two-parton process. Instead, we regard the whole system as the same temperature, though which is mass dependent.
The present work is different from the quark coalescence model [66,[76][77][78][79][80], though both the models have used the thermalization and statistics. In particular, the quark coalescence model describes classically mesonic prehadrons as quark-anti-quark clusters, and baryonic ones composed from three quarks. The present work describes both mesons and baryons as the products of two participant partons which are regarded as two energy sources.
The assumption of two participant partons discussed in the present work does not mean that the particles considered directly stem from two initial partons of the incoming nuclei. In fact, we assume the two participant partons from the violent collision system in which there is rescattering, recombination, or coalescence. The two participant partons are only regarded as two energy sources to produce a considered particle, whether it is a meson, baryon, or even a lepton [48,49]. The present work treats uniformly the production of final-state particles from the viewpoint of participant energy sources, but not the quark composition of the considered particles [66,[76][77][78][79][80].
In the two-component distribution (Equation (10)), the first component contributed by the soft excitation process is from the sea quarks. The second component contributed by the hard scattering process is from the valence quarks. This explanation is different from the Werner's picture on core-corona separation [81][82][83][84], in which core and corona are simply defined by the density of partons in a particular area of phase or coordinate-space and they distinguish between thermal and non-thermal particle production. This could also be a two-component fit based on the Tsallis function, but its relation to the two-parton dynamics pushed here is not clear. Anyhow, it is possible that the two processes can be described by a uniform method [48,49], though different functions and algorithms are used.
Although there were many papers in the past that have studied the identified particle spectra in high-energy collisions, both experimentally and phenomenologically, this work shows a new way to systemize the experimental data in AA collisions over a wide energy range from 2.7 GeV to 5.44 TeV at the parton level. We emphasize that, in this work, we have analyzed the particle p T as the root sum square of transverse momenta p t1 and p t2 of two participant partons. That is, the relation of p T = p 2 t1 + p 2 t2 is used. While, in our recent work [48,49], the relation of p T = p t1 + p t2 is used, which is considered from energy relation at mid-y for massless particle. The scenarios used in this work and our recent work are different. Based on our analyses, it is hard to judge which scenario is more reasonable.
Through the analysis of the data, we have obtained the excitation functions of some quantities, such as T and its weighted average T , T 0 and its weighted average T 0 , β T and its weighted average β T , q and its weighted average q , as well as a 0 and its weighted average a 0 . These excitation functions all show some specific laws as √ s NN increases. Although the conclusion on "onset of deconfinement" or QCD phase transition is indicated around 10 GeV or below is possibly over-interpreting the data and only using the blast-wave or Tsallis-like model is clearly not enough, the sudden change in the slope in the excitation function of T 0 is worthy of attention.

Summary and Conclusions
We summarize here our main observations and conclusions.
(a) The transverse momentum (mass) spectra of charged pions, kaons, and protons produced at mid-rapidity in central AA (Au-Au, Pb-Pb, and Xe-Xe) collisions over an energy range from 2.7 GeV to 5.44 TeV have been analyzed in this work. The experimental data measured by several collaborations are fitted satisfactorily in the framework of multisource thermal model in which the transverse momentum of identified particle is regarded as the root sum square of transverse momenta of two participant partons, where the latter obeys the revised Tsallis-like function. This treatment for the spectra of transverse momenta is novel and successful. The excitation functions of parameters such as the effective temperature, entropy index, revised index, kinetic freeze-out temperature, and transverse flow velocity are obtained. The chemical potential has obvious influence on the excitation function of kinetic freeze-out temperature at lower energy.
(b) With increasing collision energy, the entropy index increases slowly, and the revised index increases quickly and then changes slowly except for a large increase at the LHC. With increasing the particle mass, the entropy index decreases and the revised index increases obviously. The collision system discussed in this work stays approximately in the equilibrium state, and some functions based on the assumption of equilibrium can be used. The system is closer to the equilibrium state when it emits protons at lower energy, comparing with pions and kaons at higher energy. The revised index describes the bending degrees of the spectra in very low transverse momentum region. Its values for the spectra of charged pions, kaons, and protons are approximately 0.75, 1, and 1.8, respectively, at above 10 GeV and drop obviously at below 10 GeV.
(c) With increasing collision energy, the effective temperature increases clearly and monotonously, and the kinetic freeze-out temperature and transverse flow velocity increase quickly from ≈ 3 to ≈ 10 GeV and then slowly at above 10 GeV. There is a plateau from near 10 GeV to 200 GeV in the excitation functions of the latter pair. The onset energy of deconfinement phase transition from hadronic matter to QGP is connected to the special changes of excitation function of kinetic freeze-out temperature and possibly 10 GeV or slightly lower. If the plateau at the RHIC is regarded as a reflection of the formation of QGP liquid drop, the following quick increase of the excitation functions at the LHC is a reflection of higher temperature QGP liquid drop due to larger energy deposition. At kinetic freeze-out, the temperature and expansion velocity of the system increase with increasing the energy from the RHIC to LHC.

Data Availability Statement:
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.