Kinetic Freeze-Out Properties from Transverse Momentum Spectra of Pions in High Energy Proton-Proton Collisions

: Transverse momentum spectra of negative and positive pions produced at mid-(pseudo)rapidity in inelastic or non-single-diffractive proton-proton collisions over a center-of-mass energy, √ s , range from a few GeV to above 10 TeV are analyzed by the blast-wave ﬁt with Boltzmann (Tsallis) distribution. The blast-wave ﬁt results are well ﬁtting to the experimental data measured by several collaborations. In a particular superposition with Hagedorn function, both the excitation functions of kinetic freeze-out temperature ( T 0 ) of emission source and transverse ﬂow velocity ( β T ) of produced particles obtained from a given selection in the blast-wave ﬁt with Boltzmann distribution have a hill at √ s ≈ 10 GeV, a drop at dozens of GeV, and then an increase from dozens of GeV to above 10 TeV. However, both the excitation functions of T 0 and β T obtained in the blast-wave ﬁt with Tsallis distribution do not show such a complex structure, but a very low hill. In another selection for the parameters or in the superposition with the usual step function, T 0 and β T increase generally quickly from a few GeV to about 10 GeV and then slightly at above 10 GeV, there is no such the complex structure, when also studying nucleus-nucleus collisions.


Introduction
Chemical and thermal or kinetic freeze-outs are two of important stages of system evolution in high energy collisions. The excitation degrees of interacting system at the two stages are possibly different from each other. To describe different excitation degrees of interacting system at the two stages, one can use chemical and kinetic freeze-out temperatures respectively. Generally, at the stage of chemical freeze-out, the ratios of different types of particles are no longer changed, and the chemical freeze-out temperature can be obtained from the ratios of different particles in the framework of thermal model [1][2][3]. At the stage of kinetic freeze-out, the transverse momentum spectra of different particles are no longer changed, and the dissociation temperature [4] or kinetic freeze-out temperature can be obtained from the transverse momentum spectra according to the hydrodynamical model [4] and the subsequent blast-wave fit with Boltzmann distribution [5][6][7] or with Tsallis distribution [8][9][10].
It should be pointed out that the transverse momentum spectra even though in narrow range contain both the contributions of random thermal motion and transverse flow of particles. The random thermal motion and transverse flow reflect the excitation and expansion degrees of the interacting system (or emission source) respectively. To extract the kinetic freeze-out temperature from transverse 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 3.1, we summarize our main observations and conclusions.

Formalism and Method
There are two main processes of particle productions, namely the soft excitation process and the hard scattering process, in high energy collisions. For the soft excitation process, the method used in the present work is the blast-wave fit [5][6][7][8][9][10] that has wide applications in particle productions. The blast-wave fit is based on two types of distributions. One is the Boltzmann distribution [5][6][7] and another one is the Tsallis distribution [8][9][10]. As an application of the blast-wave fit, we present directly its formalisms in the following. Although the blast-wave fit has abundant connotations, we focus only our attention on the formalism of transverse momentum (p T ) distribution in which the kinetic freeze-out temperature (T 0 ) and mean transverse flow velocity (β T ) are included.
We are interested in the blast-wave fit with Boltzmann distribution in its original form. According to [5][6][7], the blast-wave fit with Boltzmann distribution results in the probability density distribution of p T to be where C 1 is the normalized constant, m T = p 2 T + m 2 0 is the transverse mass, m 0 is the rest mass, r is the radial coordinate in the thermal source, R is the maximum r which can be regarded as the transverse size of source in the case of neglecting the expansion, r/R is the relative radial position which has in fact more meanings than r and R themselves, I 0 and K 1 are the modified Bessel functions of the first and second kinds respectively, ρ = tanh −1 [β(r)] is the boost angle, β(r) = β S (r/R) n 0 is a self-similar flow profile, β S is the flow velocity on the surface, and n 0 = 2 is used in the original form [5]. There is the relation between β T and β(r). As a mean of β(r), β T = (2/R 2 ) R 0 rβ(r)dr = 2β S /(n 0 + 2). We are also interested in the blast-wave fit with Tsallis distribution in its original form. According to [8][9][10], the blast-wave fit with Tsallis distribution results in the p T distribution to be where C 2 is the normalized constant, q is an entropy index that characterizes the degree of non-equilibrium, φ denotes the azimuthal angle, and n 0 = 1 is used in the original form [8]. Because of n 0 being an insensitive quantity, the results corresponding to n 0 = 1 and 2 for the blast-wave fit with Boltzmann or Tsallis distribution are harmonious [18]. In fact, in some literature [29], n 0 is regarded as a free parameter which changes largely by several times and increases 1 in the number of free parameter, which is not our expectation in the present work. In addition, the index −1/(q − 1) used in Equation (2) can be replaced by −q/(q − 1) due to q being very close to 1. This substitution results in a small and negligible difference in the Tsallis distribution [30,31]. For a not too wide p T spectrum, the above two equations can be used to describe the p T spectrum and to extract T 0 and β T . For a wide p T spectrum, we have to consider the contribution of hard scattering process. According to the quantum chromodynamics (QCD) calculus [32][33][34], the contribution of hard scattering process is parameterized to be an inverse power-law which is the Hagedorn function [35,36], where p 0 and n are free parameters, and A is the normalization constant related to the free parameters. In literature [37][38][39][40][41][42][43], there are respectively modified Hagedorn functions and where the three normalization constants A, free parameters p 0 , and free parameters n are severally different, though the same symbols are used to avoid trivial expression. The experimental p T spectrum distributed in a wide range can be described by a superposition of the contributions of the soft excitation and hard scattering processes. If one of Equations (1) and (2) describes the contribution of the soft excitation process, one of Equations (3)-(6) describes the contribution of the hard scattering process. To describe the spectrum in a wide p T range, we can superpose a two-component superposition like this where k (1 − k) denotes the contribution fraction of the soft excitation (hard scattering) process, and f S (p T ) denotes one of Equations (1) and (2). As for the four f H (p T ), we are inclined to the first one due to its more applications. Naturally, we have the normalization condition p T max 0 f 0 (p T )dp T = 1, where p T max denotes the maximum p T . In Equation (7), the soft component contributes in low p T region and the hard component contributes in whole p T range. The two contributions overlap each other in low p T region.
According to Hagedorn's model [35], we may also use the usual step function to superpose the two functions. That is where A 1 and A 2 are constants which result in the two components to be equal to each other at p T = p 1 ≈ 2 ∼ 3 GeV/c. The contribution fraction of the soft excitation (hard excitation) process in Equation (8) f 0 (p T )dp T = 1. In Equation (8), the soft component contributes in low p T region and the hard component contributes in high p T region. The two contributions link with each other at p T = p 1 .
In some cases, the contribution of resonance production for pions in very-low p T range has to be considered. We can use a very-soft component for the p T range from 0 to 0.2∼0.3 GeV/c which covers the contribution of resonance production. Let k VS and k S denote the contribution fractions of the very-soft and soft processes respectively. Equation (7) is revised to where f VS (p T ) denotes one of Equations (1) and (2) as f S (p T ), but having smaller parameter values comparing with f S (p T ). Anyhow, both the very-soft and soft processes are belong to the soft process. Correspondingly, Equation (8) is revised to where A 1 , A 2 , and A 3 are constants which result in the two contiguous components to be equal to each other at p T = p 1 and p T = p 2 . The above two types of superpositions [Equations (7) and (8)] treat the soft and hard components by different ways in the whole p T range. Equation (7) means that the soft component contributes in a range from 0 up to 2∼3 GeV/c or a little more. The hard component contributes in the whole p T range, though the main contributor in the low p T region is the soft component and the sole contributor in the high p T region is the hard component. Equation (8) shows that the soft component contributes in a range from 0 up to p 1 , and the hard component contributes in a range from p 1 up to the maximum. The boundary of the contributions of soft and hard components is p 1 . There is no mixed region for the two components in Equation (8).
In the case of including only the soft component, Equations (7) and (8) are the same. In the case of including both the soft and hard components, their common parameters such as T 0 , β T , p 0 , and n should be severally to have small differences from each other. To avoid large differences, we should select the experimental data in a narrow p T range. In addition, the very-soft component in Equations (9) and (10) does not need to consider in some cases due to the fact that the spectrum in very-low p T range are possibly not available in experiment. Then, Equations (9) and (10) are degenerated to Equations (7) and (8) respectively in some cases.
In the fit process, firstly, we use Equation (7) to extract the related parameters, where f S (p T ) and f H (p T ) are exactly Equations (1) or (2) and (3) respectively. Regardless of Equation (1) or (2) is regarded as f S (p T ), the situation is similar due to Equations (1) and (2) being harmonious in trends of parameters [18], though one more parameter (the entropy index q) is needed in Equation (2). Secondly, we use Equation (8) to extract the related parameters as comparisons with those from Equations (7). In the case of Equations (7) and (8) being not suitable, Equations (9) and (10) can be used due to the very-low component from resonance decays being included in the first items of the two fits. Meanwhile, the contribution of resonance decays in the low component is naturally included in the second items of Equations (9) and (10) or in the first items of Equations (7) and (8). Thus, the contribution of resonance decays which is available in the very-low or low component in experiment is naturally considered by us.
It should be pointed out that although Equations (9) and (10) are not used in the final fits in the present work, we hope to keep them here due to the fact that we need them to give further explanations. From theses explanations, one can see the relations between Equations (9)/(10) and (7)/(8), as well as possible applications of Equations (9) and (10). Figure 1 shows the transverse momentum spectra of π − and π + produced at mid-(pseudo)rapidity in pp collisions at high center-of-mass energies, where different mid-(pseudo)rapidity (y or η) intervals and energies ( √ s) are marked in the panels. Different forms of the spectra are used due to different Collaborations, where N, E, p, σ, and N EV denote the particle number, energy, momentum, cross-section, and event number, respectively. The closed and open symbols presented in panels (a)-(e) represent the data of π − and π + measured by the NA61/SHINE [24], PHENIX [25], STAR [6] [44][45][46][47][48] given no data on the p T spectra of π − and π + . These data show wider p T spectra of charged particles to be used in our foreseen studies. In some cases, different amounts marked in the panels are used to scale the data for clarity. We have fitted the data by two sets of parameter values in Equation (7) so that we can see the fluctuations of parameter values. The blue solid (dotted) curves are our results for π − (π + ) spectra fitted by Equation (7) through Equations (1) and (3), and the blue dashed (dot-dashed) curves are our results for π − (π + ) spectra fitted by Equation (7) through Equations (2) and (3), by the first set of parameter values, in which k is taken to close to 0.5 as much as possible. The black curves are our results fitted by the second set of parameter values for comparison, in which k is taken to close to 1 as much as possible. The values of free parameters (T 0 , β T , q if available, k, p 0 , and n), normalization constant (N 0 ), χ 2 , and number of degrees of freedom (ndof) corresponding to the curves in Figure 1 are listed in Tables 1 and 2, where the errors of fit parameters are obtained by the statistical simulation method [49], no matter what χ 2 /ndof is. The parameter values presented in terms of value 1 /value 2 denote respectively the first and second sets of parameter values in Equation (7) through Equations (1)-(3) in which k = 1. One can see that Equation (7) with two sets of parameter values describes the p T spectra at mid-(pseudo)rapidity in pp collisions over an energy range from a few GeV to above 10 TeV. The blast-wave fit with Boltzmann distribution and with Tsallis distribution presents similar results. The free parameters show some laws in the considered energy range. For a given parameter, its fluctuation at given energy is obvious in some cases. Because of the data being not available in very-low p T range, Equation (9) is not used in the fit.

Comparison with Data by Equation (7)
It should be noted that, from the fit process we know that, the Tsallis expression, Equation (2), having a polynomial behavior at large p T could be a better description for the spectra at large values of p T , where the Boltzmann expression, Equation (1) does not work possibly at large p T if we use the same parameters T 0 and β T . In some cases, for the spectra in a wide p T range, a single Tsallis expression is suitable, and two-or three-component Boltzmann expression is needed. Our previous work [50] studied both the Tsallis and Boltzmann distributions without flow effect and also confirms this issue. Indeed, the Tsallis description is better than the Boltzmann one, though the former has one more parameter q. In fact, the introduction of the entropy index q in the Tsallis description has the meaning of reality. As discussed in Section 2, q characterizes the degree of non-equilibrium. In addition, when q → 1, the Tsallis description degenerates to the Boltzmann one.  (e) Figure 1. Transverse momentum spectra of π − and π + produced at mid-(pseudo)rapidity in pp collisions at high energies, where the mid-(pseudo)rapidity intervals and energies are marked in the panels. The symbols presented in panels (a-e) represent the data of NA61/SHINE [24], PHENIX [25], STAR [6] In some cases, different amounts marked in the panels are used to scale the data for clarity. The blue solid (dotted) curves are our results for π − (π + ) spectra fitted by Equation (7) through Equations (1) and (3), and the blue dashed (dot-dashed) curves are our results for π − (π + ) spectra fitted by Equation (7) through Equations (2) and (3), by the first set of parameter values. The results by the second set of parameters (if available) are presented by the black curves. Table 1. Values of free parameters (T 0 , β T , k, p 0 , and n), normalization constant (N 0 ), χ 2 , and ndof corresponding to the solid (dotted) curves for π − (π + ) spectra in Figure 1 in which different data are measured in different mid-(pseudo)rapidity intervals at different energies by different Collaborations. The values presented in terms of value 1 /value 2 denote respectively the first and second sets of parameter values in Equation (7) through Equations (1) and (3) in which k = 1.  Table 2. Values of free parameters (T 0 , β T , q k, p 0 , and n), normalization constant (N 0 ), χ 2 , and ndof corresponding to the dashed (dot-dashed) curves for π − (π + ) spectra in Figure 1 in which different data are measured in different mid-(pseudo)rapidity intervals at different energies by different Collaborations. The values presented in terms of value 1 /value 2 denote respectively the first and second sets of parameter values in Equation (7) through Equations (2) and (3) in which k = 1. To see clearly the excitation functions of free parameters, Figure 2a-e show the dependences of T 0 , β T , p 0 , n, and k on √ s, respectively. The blue and black closed and open symbols represent the parameter values corresponding to π − and π + respectively, which are listed in Tables 1 and 2. The blue circles (black squares) represent the first set of parameter values obtained from Equation (7) through Equations (1)-(3). The blue asterisks (black triangles) represent the second set of parameter values obtained from Equation (7) through Equations (1)-(3). One can see that the difference between the results of π − and π + is not obvious. In the excitation functions of the first set of T 0 and β T obtained from the blast-wave fit with Boltzmann distribution, there are a hill at √ s ≈ 10 GeV, a drop at dozens of GeV, and then an increase from dozens of GeV to above 10 TeV. In the excitation functions of the first set of T 0 and β T obtained from the blast-wave fit with Tsallis distribution, there is no the complex structure, but a very low hill. In the excitation functions of the second set of T 0 and β T , there is a slight increase from about 10 GeV to above 10 TeV. In Equation (7) contained the blast-wave fit with both distributions, in the excitation functions of p 0 and n, there are a slight decrease and increase respectively in the case of the hard component being available. The excitation function of k shows that the contribution (1 − k) of hard component slightly increases from dozens of GeV to above 10 TeV, and it has no contribution at around 10 GeV. At given energies, the fluctuations in a given parameter result in different excitation functions due to different selections. As a comparison, the red asterisks (green triangles) in Figure 2 represent the results from AA collisions which are discussed in detail in the Appendix A. One can see that the results from AA collisions approach to those from pp collisions with the second set of parameters, though only the soft component is used in most cases.
Indeed, √ s NN ≈ 10 GeV is a special energy for AA collisions as indicated by Cleymans [51].
The present work shows that √ s ≈ 10 GeV is also a special energy for pp collisions. In particular, there is a hill in the excitation functions of T 0 and β T in pp collisions due to a given selection of the parameters. At this energy (11 GeV more specifically [51]), the final state has the highest net baryon density, a transition from a baryon-dominated to a meson-dominated final state takes place, and the ratios of strange particles to mesons show clear and pronounced maximums [51]. These properties result in this special energy.
At 11 GeV, the chemical freeze-out temperature in AA collisions is about 151 MeV [51], and the present work shows that the kinetic freeze-out temperature in pp collisions is about 105 MeV, extracted from the blast-wave fit with Boltzmann distribution. If we do not consider the difference between AA and pp collisions, though cold nuclear effect exists in AA collisions, the chemical freeze-out happens obviously earlier than the kinetic one. According to an ideal fluid consideration, the time evolution of temperature follows T f = T i (τ i /τ f ) 1/3 , where T i (= 300 MeV) and τ i (= 1 fm) are the initial temperature and proper time respectively [52,53], and T f and τ f denote the final temperature and time respectively, the chemical and kinetic freeze-outs happen at 7.8 and 23.3 fm respectively. It should be noted that in the calculation of τ f , the Lorentz factor is not considered. If we consider the mean Lorentz factor (γ ≈ 5-6) of charged pions in the rest frame of emission source [15][16][17][18], the value of freeze-out time will be smaller.  Tables 1 and 2.
The blue circles (black squares) represent the first set of parameter values obtained from Equation (7) through Equations (1)-(3). The blue asterisks (black triangles) represent the second set of parameter values obtained from Equation (7) through Equations (1)-(3). The red asterisks (green triangles) represent the parameter values from AA collisions for comparisons, which are listed in Tables A1 and A2 in the Appendix A.
Strictly, T 0 (β T ) obtained from the pion spectra in the present work is less than that averaged by weighting the yields of pions, kaons, protons, and other light particles. Fortunately, the fraction of the pion yield in high energy collisions are major (≈ 85%). The parameters and their tendencies obtained from the pion spectra are similar to those obtained from the spectra of all light particles. To study the excitation functions of T 0 and β T , it does not matter if we use the spectra of pions instead of all light particles.
It should be noted that the main parameters T 0 and β T are correlated in some way. Although the excitation functions of T 0 (β T ) which are acceptable in the fit process are not sole, their tendencies are harmonious in most cases, in particular in the energy range from the RHIC to LHC. Combining with our previous work [18], we could say that there is a slight (≈ 10%) increase in the excitation function of T 0 and an obvious (≈ 35%) increase in the excitation function of β T from the RHIC to LHC. At least, the excitation functions of T 0 and β T do not decrease from the RHIC to LHC. However, the excitation function of T 0 from low to high energies is not always incremental or invariant, though the excitation function of β T has the trend of increase in general. For example, In [4], T 0 slowly decreases as √ s increases from 23 GeV to 1.8 TeV, and β T slowly increases with √ s. In [19,20], T 0 has no obvious change and β T has a slight (≈ 10%) increase from the RHIC to LHC. In [21], T 0 has a slight (≈ 9%) increase and β T has a large (≈ 65%) increase from the RHIC to LHC. In [22,23], T 0 has a slight (≈ 5%) decrease from the RHIC to LHC and β T increases by ≈ 20% from 39 to 200 GeV. It is convinced that β T increases from the RHIC to LHC, though the situation of T 0 is doubtful.
Although some works [54][55][56][57] reported a decrease of T 0 and an increase of β T from the RHIC to LHC, our re-scans on their plots show a different situation of T 0 . For example, in [54], our re-scans show that T 0 has no obvious change and β T has a slight (≈ 9%) increase from the top RHIC to LHC, though there is an obvious hill or there is an increase by ≈ 30% in T 0 in 5-40 GeV comparting with that at the RHIC. Authors in [55] shows similar results to [54] with the almost invariant T 0 from the top RHIC to LHC, an increase by ≈ 28% in T 0 in 7-40 GeV comparing with that at the top RHIC, and an increase by ≈ 8% in β T comparing with that at the top RHIC. Authors in [56,57] shows similar result to [54,55] on T 0 , though the excitation function of β T is not available.
In some cases, the correlation between T 0 and β T are not negative, though some works [54,55] show negative correlation over a wide energy range. For a give p T spectrum, it seems that a larger T 0 corresponds to a smaller β T , which shows a negative correlation. However, this negative correlation is not sole case. In fact, a couple of suitable T 0 and β T can fit a given p T spectrum. A series of p T spectra at different energies possibly show a positive correlation between T 0 and β T , or independent of T 0 on β T , in a narrow energy range. Very recently, authors in [58] shows approximately independent of T 0 on β T in pp collisions at √ s = 7 TeV, and negative correlations in p-Pb collisions at √ s NN = 5.02 TeV and in Pb-Pb collisions at √ s NN = 2.76 TeV, for different average charged-particle multiplicity densities, i.e., for different centrality classes. These results are partly in agreement with our results. It seems that the correlation between T 0 and β T is an open question at present, though some researchers think that there is a negative correlation between T 0 and β T . In our opinion, the type of correlation between T 0 and β T depends on three factors, that is the choices of fitted region in low and medium p T range, fixed or changeable n 0 , sensitivity of β T on centrality. Indeed, more studies are needed in the near future.
To study further the behaviors of parameters, Figure 3 shows the excitation functions of (a)(b) mean p T ( p T ) and (c)(d) ratio of root-mean-square p T ( p 2 T ) to with χ 2 /ndof = 54/82. From the line one can see that the behavior of p T . In particular, with the increase of ln √ s and including the contribution of second component, p T increases approximately linearly. As a comparison, the red asterisks (green triangles) in Figure 3 represent the results for AA collisions, which are indirectly obtained according to the parameter values listed in Tables A1 and A2. One can see that the results for AA collisions are greater than those for pp collisions at around 10 GeV and above.
The quantities p T and T i are very important to understand the excitation degree of interacting system. As for the right panel in Figure 3 which is for the two-component, the incremental trend for p T and T i with the increase of √ s ( √ s NN ) is a natural result due to more energy deposition at higher energy. Although p T and T i are obtained from the parameter values listed in Tables 1 and 2 (A1 and A2), they are independent of fits or models. More investigations on the excitation functions of p T and T i are needed due to their importance.  (8) To discuss further, for comparisons with the results from Equation (7), we reanalyze the spectra by Equation (8) and study the trends of new parameters. Figure 4 is the same as Figure 1, but showing the results fitted by Equation (8) through Equations (1) and (3) and through Equations (2) and (3) respectively. For Equation (8) Table 3. The related parameters are shown in Figure 5 and the leading-out parameters are shown in Figure 6, which are the same as Figures 2 and 3 respectively, and with only one set of parameter values for Equation (8) through Equations (1) and (3). In particular, k in Figure 5 is obtained by k = p 1 0 A 1 f S (p T )dp T due to f 0 (p T ) is normalized to 1, as discussed following Equation (8). The lines in Figure 6a

Comparison with Data by Equation
and with χ 2 /ndof = 55/61 and 28/61 respectively, though the linear relationships between the parameters and ln √ s may be not the best fitting functions. Similar to Figures 2 and 3, the results for AA collisions from Equation (8) are also presented in Figures 5 and 6 for comparisons. One can see the similarity in both pp and AA collisions in the considered energy range.
From Figure 4 one can see that Equation (8) fits similarly good the data as Equation (7). Because of the data being not available in very-low p T range, Equation (10) is not used in the fit. T 0 and β T in Figure 5 increase slightly from a few GeV to above 10 TeV with some fluctuations in some cases, which is partly similar to those in Figure 2. Other parameters in Figure 5 show somehow similar trends toFigure 2 with some differences. The left panels in Figures 3 and 6 Table 3. Values of T 0 , β T , k, p 0 , n, N 0 , χ 2 , and ndof corresponding to the solid (dotted) curves for π − (π + ) spectra in Figure 4, where Equation (8) Table 4. Values of T 0 , β T , q, k, p 0 , n, N 0 , χ 2 , and ndof corresponding to the dashed (dot-dashed) curves for π − (π + ) spectra in Figure 4, where Equation (8) through Equations (2) and (3) is used.    (2) and (3) with two sets of parameter values. The blue solid (dotted) curves are the results for π − (π + ) spectra fitted by Equation (8) through Equations (1) and (3), and the blue and black dashed (dot-dashed) curves are the results for π − (π + ) spectra fitted by Equation (8) through Equations (2) and (3). Combining with our recent work [62], it is shown the similarity in pp and AA collisions, though AA collisions appear larger T 0 , β T , p T , and T i in most cases. Moreover, it is well seen that Tsallis does not distinguish well between the data in pp and AA collisions [63][64][65]. Indeed, at high energy, both pp and AA collisions produce many particles and obey statistical law. In addition, pp collisions are the basic sub-process in AA collisions. It is natural that pp and AA collisions show similar law. However, concerning around 10-20 GeV change, this is expected as soon as QCD effects/calculations may not be directly applicable below this point, plus seem to be smoothed away by other processes in AA collisions. This is well visible as a clear difference as shown in [63][64][65], where the data in pp collisions goes well with the data in electron-positron collisions, while the data in AA collisions are different.
The differences between Equations (7) and (8) are obviously, though the similar components are used in them. In our recent works [18,66], Equations (7) and (8) are used respectively. Although there is correlation in the extraction of parameters, a smooth curve can be easily obtained by Equation (7). Although it is not easy to obtain a smooth curve at the point of junction, there is no or less correlation in the extraction of parameters by Equation (8). In consideration of obtaining a set of parameters with least correlation, we are inclined to use Equation (8) to extract the related parameters. This inclining results in Equation (8) to separate determinedly the contributions of soft and hard processes.
It should be noted that the system of pp collisions at low energy is probably not in thermal equilibrium or local thermal equilibriums due to low multiplicity, which is not the case of the present work. In fact, the present work treats pp collisions at high energies in which the multiplicities in most cases are not too low. In addition, related review work [67] shows that small system also appears similar collective behavior to AA collisions. This renders that the idea of local thermal equilibrium is suitable to high energy pp collisions for which the blast-wave fit can be used, though QGP is not expected to form in minimum-bias events.
Although we have used the blast-wave fit in the superposition function with two components in which the second component is an inverse power-law, the blast-wave part is not necessary for fitting process itself. In fact, the superposition of (two-)Boltzmann (or Tsallis) distribution and inverse power-law can fit the data in most cases [30,31,65]. In particular, in our very recent work [68], we used the Tsallis-Pareto-type function [28, 69,70] to fit p T spectra in a wide range, in which there is no boosted part. The merits of the boosted part are that some additional quantities such as T 0 and β T can be obtained and physics picture is more abundant.
To avoid the dependences of T 0 and β T on fits or models, one can use p T to describe synchronously T 0 and β T . Generally, p T is independent of fits or models, though it can be calculated from fits or models. Averagely, the contribution of one participant in each binary collisions is p T /2 which is resulted from both the thermal motion and flow effect. Let k 0 denote the contribution fraction of thermal motion. The contribution fraction of flow effect is naturally 1 − k 0 . Thus, we define T 0 = k 0 p T /2 and β T = (1 − k 0 ) p T /2m 0 γ. As a free parameter, k 0 depends on collision energy, which is needed to study further.

More Discussions
Before summary and conclusions, we would like to underline the preponderance of the present work. Comparing with PYTHIA or other perturbative QCD simulation tools [71][72][73][74], the present work is simpler and more applicative in obtaining the excitation functions of T 0 and β T , though the usual fitting method is used. In a recent work [75], it was pointed out that the PYTHIA Monte Carlo [71,72] disagrees with some data, and the two-component (soft+hard) model describes the data accurately and comprehensively, though the two-component model used in [75] is different from the fit used in the present work. In addition, the present work is a data-driven reanalysis based on some physics considerations, but not a simple fit to the data. From the data-driven reanalysis, the excitation functions of some quantities have been obtained.
From the excitation functions, one can see some complex structures which are useful to study the properties of particle production and system evolution at different energies. In particular, in the energy range around 10 GeV, the excitation functions have transition which implies the phase of interaction matter had changed. In addition, by using the two-component fit, the present work also presents a new method to extract the contribution fraction of hard component. One can see that this contribution fraction is 0 in the energy range around 10 GeV. This implies that the interactions in the energy range around 10 GeV have only soft component, and that above 10 GeV have both soft and hard components. The interaction mechanisms in the two energy ranges are different.
We would like to emphasize that the aim of the present work is to look for possible signatures of a transition from baryon-dominated to meson-dominated hadron production mechanism by fitting pion p T spectra at mid-(pseudo)rapidity in pp collisions which also show collective phenomenon [67]. The main conclusion regards a possible indication of such an effect at about 10-20 GeV visible in a drop of the temperature extracted using a blast-wave model with Boltzmann distribution to model the soft excitation mode. Such a conclusion seems not straightforward since the results are strongly biased by different p T max used to fit data at different energies. In fact, the results are less affected by p T max due to the fact that the parameters are mainly determined by the spectra in low and intermediate p T regions. Larger p T max does not change the trend of inverse power-law, which does not affect obviously the parameters. Although p T max has no obvious influence on the parameters, we have used p T max = 5 Gev/c in our calculation.
It should be noted that the blast-wave analysis is known to be sensitive to the selected p T range in the case of using a not too wide and local one such as 1-2 GeV/c. To avoid this dependence, we have used a wide enough p T range from 0 to a large value, but not a narrow and local one. In addition, looking at the fit results when using a Tsallis function to describe soft excitation processes such a drop in the k-parameter is no longer visible which seems to mean that the fit is less sensitive to the hard component and also less sensitive to p T max . Using the Tsallis function assumption also the drop in the temperature is no longer visible which means that the effect on the temperature seen with the Boltzmann assumption seems to be model dependent. In fact, the parameters discussed in the present work are indeed model dependent. We hope to extract model independent parameters in the near future. Our definition T 0 = k 0 p T /2 and β T = (1 − k 0 ) p T /2m 0 γ are a possible choice for the model independent parameters.
On the other hand, although comparing Boltzmann with Tsallis is an interesting physical problem, which could not be understood without another comparison with the generic super-statistics, as the one implemented in particle productions [76][77][78][79]. The latter-in contrary to Boltzmann and Tsallis-lets the system alone judge about its statistical nature, whether extensive or non-extensive (equilibrium or non-equilibrium). As pointed out in [76][77][78][79], the particle production at BES energies is likely a non-extensive process but not necessarily Boltzmann or Tsallis type. Indeed, further study on particle production in high energy collisions is needed in the future.
In particular, larger T 0 means higher excitation degree and shorter lifetime of the fireball formed in high energy collisions. At the LHC energy, the fireball should have higher excitation degree comparing with that at the top RHIC energy, though longer lifetime is possible at the LHC energy [55,80]. As a result of competition between excitation degree and lifetime, T 0 shows the trend of increase with energy in the present work. Meanwhile, our result on larger β T means quicker expansion at the LHC energy. These trends are harmonious with those of p T and T i as well as other method such as the alternative method [16,18,66] by which T 0 and β T are also obtained.

Summary and Conclusions
The transverse momentum spectra of π − and π + produced at mid-(pseudo)rapidity in pp collisions over an energy range from a few GeV to above 10 TeV have been analyzed by the superposition of the blast-wave fit with Boltzmann distribution or with Tsallis distribution and the inverse power-law (Hagedorn function). The fit results are well fitting to the experimental data of NA61/SHINE, PHENIX, STAR, ALICE, and CMS Collaborations. The values of related parameters are extracted from the fit process and the excitation functions of parameters are obtained.
In the particular superposition Equation (7) and with a given selection, both excitation functions of T 0 and β T obtained from the blast-wave fit with Boltzmann distribution show a hill at √ s ≈ 10 GeV, a drop at dozens of GeV, and an increase from dozens of GeV to above 10 TeV. The mentioned two excitation functions obtained from the blast-wave fit with Tsallis distribution does not show such a complex structure, but a very low hill. In another selection for the parameters in Equation (7) or in the superposition Equation (8), T 0 and β T increase generally quickly from a few GeV to about 10 GeV and then slightly at above 10 GeV. There is no the complex structure, too. In both superpositions, the excitation function of p 0 (n) shows a slight decrease (increase) in the case of the hard component being available. From the RHIC to LHC, there is a positive (negative) correlation between T 0 and β T (p 0 and n). The contribution of hard component slightly increases from dozens of GeV to above 10 TeV, and it has no contribution at around 10 GeV.
In the case of considering the two components together, the mean transverse momentum and the initial temperature increase obviously with the increase of logarithmic collision energy in the considered energy range. From a few GeV to above 10 TeV, the collision system takes place possibly a transition at around 10 GeV, where the transition from a baryon-dominated to a meson-dominated final state takes place. No matter what a structure appears, the energy range around 10 GeV is a special one due to the slope of T 0 excitation function having large variation. Indeed, the mentioned energy range is needed further study in the future. Acknowledgments: Communications from Edward K. Sarkisyan-Grinbaum are highly acknowledged.

Conflicts of Interest:
The authors declare that there are no conflicts of interest regarding the publication of this paper. 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.

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.
Compliance with Ethical Standards: The authors declare that they are in compliance with ethical standards regarding the content of this paper.

Appendix A. Fit Results from AA Collisions
Although our previous work [62] studied the excitation functions of T 0 and β T in AA collisions, different fit functions were used there. To make a comparison with the results from pp collisions in the present work, we have to use the same fit functions, i.e., Equations (7) and (8), to re-fit the spectra from AA collisions. Figure A1 shows the spectra of π − and π + produced in mid-rapidity range in central Cu-Cu, Au-Au, and Pb-Pb collisions at high energies. Panels (a)-(f) represent the data by various symbols measured by the FOPI [81], STAR [55], STAR [7], PHENIX [82], STAR [6], and ALICE [20] Collaborations, respectively, where for 2.24 GeV Au-Au collisions in panel (a) only the spectrum of π − is available. As the same as Figure 1, the solid (dotted) curves are our results for π − (π + ) spectra fitted by Equation (7) through Equations (1) and (3), and the dashed (dot-dashed) curves are our results for π − (π + ) spectra fitted by Equation (7) through Equations (2) and (3), where only one set of parameters is used. The values of parameters are listed in Tables A1 and A2 with χ 2 and ndof, where the values of parameters are used directly in Figure 2 and indirectly in Figure 3. The fit results based on Equation (7) are well fitting to the experimental data measured in central AA collisions by the international collaborations [6,7,20,55,81,82]. Figure A2 is the same as Figure A1, but it shows the fit results from Equation (8) through Equations (1) and (3) as well as through Equations (2) and (3). As the same as Figure 4, the solid (dotted) curves are the results for π − (π + ) spectra from Equation (8) through Equations (1) and (3), and the dashed (dot-dashed) curves are the results for π − (π + ) spectra from Equation (8) through Equations (2) and (3). The values of parameters are listed in Tables A3 and A4 with χ 2 and ndof, where the values of parameters are used directly in Figure 5 and indirectly in Figure 6. The fit results based on Equation (8) are well fitting to the experimental data measured in central AA collisions by the international collaborations [6,7,20,55,81,82]. It should be noted that Tables A3 and A4 are nearly the  same as Tables A1 and A2 respectively, if not equal in error. The reason is that narrow p T ranges are used, in which the first component plays complete or main rule. The difference between Equations (7) and (8) is then not obvious.  Figure A1. Same as Figure 1, but showing the results from central AA collisions. Panels (a-f) represent the data by various symbols measured by the FOPI [81], STAR [55], STAR [7], PHENIX [82], STAR [6], and ALICE [20] Collaborations, respectively.  Table A1. Values of T 0 , β T , k, p 0 , n, N 0 , χ 2 , and ndof corresponding to the solid (dotted) curves for π − (π + ) spectra in Figure A1.  Table A2. Values of T 0 , β T , q, k, p 0 , n, N 0 , χ 2 , and ndof corresponding to the dashed (dot-dashed) curves for π − (π + ) spectra in Figure A1.  Table A3. Values of T 0 , β T , k, p 0 , n, N 0 , χ 2 , and ndof corresponding to the solid (dotted) curves for π − (π + ) spectra in Figure A2.  Table A4. Values of T 0 , β T , q, k, p 0 , n, N 0 , χ 2 , and ndof corresponding to the dashed (dot-dashed) curves for π − (π + ) spectra in Figure A2.