The universal scaling of kinetic freeze-out parameters across different collision systems at the LHC energy

In this paper, we perform the Tsallis Blast-Wave analysis on the transverse momentum spectra of identified hadrons produced in a wide range of collision systems at the Large Hadron Collider (LHC) including pp, pPb, XeXe and PbPb collisions. The kinetic freeze-out properties are investigated across these systems varying with the event multiplicity. We find that the extracted kinetic freeze-out temperature, radial flow velocity and the non-extensive parameter exhibit a universal scaling behavior for these systems with very different geometric size, especially when the independent baryon Tsallis non-extensive parameter is considered. This universality may indicate the existence of a unified partonic evolution stage in different collision systems at the LHC energies.

In this paper, we perform the Tsallis Blast-Wave analysis on the transverse momentum spectra of identified hadrons produced in a wide range of collision systems at the Large Hadron Collider (LHC) including pp, pPb, XeXe and PbPb collisions. The kinetic freeze-out properties are investigated across these systems varying with the event multiplicity. We find that the extracted kinetic freezeout temperature, radial flow velocity and the non-extensive parameter exhibit a universal scaling behavior for these systems with very different geometric size, especially when the independent baryon Tsallis non-extensive parameter is considered. This universality may indicate the existence of a unified partonic evolution stage in different collision systems at the LHC energies.

I. INTRODUCTION
It is indicated by lattice quantum chromodynamics (QCD) calculations that a deconfined quark-gluon plasma (QGP) state of nuclear matter might exist with high enough temperatures and densities [1]. These extreme conditions are believed to be achieved in relativistic nucleus nucleus (AA) collisions, in which the QCD parton degrees of freedom are released from the nucleons and interact with each other creating the thermalized QGP medium. The deconfined parton matter expands rapidly under the thermal pressure against the surrounding vacuum until the temperature of the system drops below which the partons are converted to hadrons, developing substantial expansion flow in the parton evolution stage [2]. After further rescatterings between the produced hadronic objects for a while, the system becomes so dilute that all particle interactions are ceased [3]. The evolving information of the medium usually described by the relativistic fluid hydrodynamics [4,5] is thus encoded in the momentum distributions of the final state particles. The transverse momentum (p T ) spectra of identified hadrons can be utilized to extract the system properties even in the early stage of the evolution within the Blast-Wave global analysis framework [6,7]. In the Blast-Wave model, particle spectra are obtained within a hydrodynamic framework through the thermal particle emissions from a freeze-out surface of the flowing medium with the kinetic freeze-out temperature T and the radial expansion flow velocity profile β(r).
The Boltzmann-Gibbs Blast-Wave (BGBW) model has been widely used to study the p T distributions of identified hadrons in heavy ion collisions. This approach assumes local equilibrium of the system so that a Boltzmann distribution can be applied to describe the emitting particles in the local rest frame of the expanding fluid. Considering the equilibrium assumption generally fails at * zhengliang@cug.edu.cn high p T , this BGBW framework is believed to work only in the low p T region and sensitive to the choice of the fitted p T range [8,9]. The Blast-Wave model is further developed with the inclusion of the Tsallis statistics [10], by assuming the Tsallis distributions for the emitting particles rather than the exponential type distributions used in the BGBW model. The Tsallis Blast-Wave (TBW) model is expected to account for the non-equilibrium or hard scattering effects by incorporating the non-extensive parameter q and describe the final state hadron spectra for an extended p T range [11][12][13][14][15][16]. Being successful in understanding a wide range of complex systems and utilized to study various high energy collisions [17][18][19], the physical interpretation of the non-extensive parameter q in Tsallis statistical method is still under discussion. It is found that the hard scattering process in high energy collisions can deliver a power law distribution of high p T hadrons from jet fragmentations with the power index related to the non-extensive parameter q [20]. This q parameter characterizes the degree of deviation from the equilibrium assumption. If q approaches unity, the TBW model returns to the BGBW model.
The energy dependence of the radial flow and the kinetic freeze-out temperature has been systematically studied [19,[21][22][23][24][25] With the unique capabilities to account for the non-equilibrium effects, the TBW model is also believed to be useful in describing the hadron productions in small systems like proton proton (pp), proton nucleus (pA) and peripheral AA collisions [26][27][28]. The recent collectivity like behaviors observed in small collision systems receives a lot of discussions on whether the QGP matter is formed in these collisions [29][30][31][32][33][34][35]. The universal strangeness enhancement effect scaled with the event multiplicity found at the LHC energy suggests a unified mechanism might induce the enhanced multi-strange hadron productions in both small and large systems [36][37][38]. Searching for the universality in other collectivity related effects like the radial expansion flow velocity and the freeze-out temperature across different collision systems and studying its system size dependence arXiv:2211.08103v2 [hep-ph] 6 May 2023 would be important to understand the appearance of the QGP like effects in small systems. It is speculated that initial energy density induced hot spots caused by Color Glass Condensate formalism may generate strong collective flow and sizable temperature fluctuations especially in small systems. The imprints of the initial fluctuations survived in the final state hadron momentum spectra at low and intermediate p T region are expected to be captured in the TBW model with the non-extensive parameter q. The temperature and flow velocity extracted from the TBW model are correlated with the q parameter by the shear and bulk viscosity in linear and quadratic forms [12,39]. Investigating the kinetic freeze-out features dependent on the event multiplicity may shed a light on the understanding of the origin of the collectivity like behavior observed in small systems.
The extracted parameters T , β and q will be systematically compared across different systems with and without considering the separate non-equilibrium parameter q for baryons. The rest of this paper is organized as follows: we illustrate the implementation of the Tsallis Blast-Wave model approach in Sec. II. The results of the extracted kinetic freeze-out parameters are presented in Sec. III. The major conclusions and discussions are summarized in Sec. IV.

II. THE TSALLIS BLAST-WAVE MODEL
The BGBW model derived from the hydrodynamic framework has been widely used to fit the transverse momentum distribution data in heavy ion collisions. It is well known that the BGBW model follows the local thermal equilibrium assumption and only works in the low p T region of the momentum spectra. The Tsallis distribution smoothly connects the power law type jet induced high p T part and the exponential type hydrodynamics dominant low p T region of the momentum distributions. With the inclusion of the non-extensive parameter, the TBW model is expected to describe the evolution from pp collisions to the central AA collisions in a consistent way. It is straightforward to implement the TBW model by replacing the particle emission sources in the local rest frame of the fluid cell from the Boltzmann distribution to the Tsallis distribution.
The invariant differential particle yield for a hadron with mass m in TBW model can be written in the form of T is the global temperature of the expanding thermal sources from which particles are emitted. R is the hard sphere edge along the transverse radial direction. A denotes the normalization constant. m T = p 2 T + m 2 is the transverse mass of a particle. y s represents the source rapidity. y b gives the beam rapidity. φ is the particle emission angle with respect to the flow velocity. ρ = tanh −1 β(r) is the radial flow profile obtained in a self-similar way with the transverse flow velocity parametrized as β(r) = β S ( r R ) n under the longitudinal boost invariant assumption. β(r) is defined in the range 0 r R with the surface velocity β S at the edge of the fireball hard sphere R and the flow profile index n. The average transverse flow velocity can be expressed as β = β S · 2/(2 + n). The relevant model parameters can be extracted based on Eq. 1 using the least χ 2 method.
In this study, the value of n is fixed to 1 for the consideration of a linear velocity profile, thus one gets β = 2/3β S . It has been argued that in small collision systems, baryon number might play an important role in hadron production and the characteristic grouping of meson and baryon sector are needed to provide better description to the p T spectra in TBW fits [26,27]. In the rest of this paper, we refer to the fits including all of the mesons and baryons with a combined non-extensive parameter q as the default TBW fit and those with four fit parameters using different q for meson (q M ) and baryon (q B ) as the TBW4 fit. In this paper, we perform both the TBW and TBW4 fits to the charged pion, kaon, and proton p T spectra at the LHC energies measured by the ALICE collaboration to extract the kinetic freeze-out parameters.

A. Transverse momentum spectra
This section compares TBW and TBW4 model fits of the transverse momentum spectra for charged pion, kaon and proton in different collision systems at the LHC energies. We restrict ourselves to the comparable spectra range p T < 3 GeV/c across all the collision systems to extract the bulk features. The average flow velocity β is required to be less than 2/3 and greater than 0 during the fit procedures to achieve a better convergence and get rid of the non-physical parameter space. For some of the very peripheral collisions, the flow parameters are very close to the boundary. We fix the flow velocity parameter β = 0 to reduce the uncertainty of the fits for those peripheral bins. It has been verified that fixing the flow velocity to zero in these peripheral bins does not change the values of the other extracted model parameters.
We present the transverse momentum spectra fits based on the default TBW model analysis in Fig. 1. Only examples for four collision systems from central to peripheral centralities are demonstrated in this figure, more details about the extracted fitting parameter values can be found in Tab. I. The black solid circles, blue solid squares and red open circles represent the experimental data in each panel. The solid lines mark the corresponding fit functions to each particle species. The results of the extracted model parameters and the χ 2 /nDoF values are also displayed in the figure. It is shown that a stronger radial flow is expected in central collisions compared to that in peripheral collisions for all systems. The freeze-out temperature mildly changes with the centrality in each collision system. The non-extensive parameter grows from peripheral to central collisions in pp collisions but drops very fast in AA collisions approaching 1 when the system is reaching equilibrium. The χ 2 /nDoF values become quite large in both very small systems like pp collisions at √ s = 7 TeV and very large systems like PbPb collisions at √ s N N = 5.02 TeV. The deviations of the fits to experimental data divided by experimental uncertainties (usually defined as pull = (f it − data)/(data error)) are shown in Fig. 2. The pull distributions can be used to quantify the agreement between the model fits and the experimental data. In all systems, the fit results seem to be consistent with data in the intermediate p T region around 1.5 GeV/c, but sizable deviations can be found generally in both smaller and higher p T regime. It is speculated that the deviations at low p T part may come from contributions due to the resonance decay effects. For a comparison, we also show the transverse momentum spectra and pull distributions based on the TBW4 fits in Fig. 3 and Fig. 4, more details about the extracted fitting parameters can be found in Tab. II. Modifying the TBW model to enable the independent q B parameter for baryons has reached high quality fits for small systems, while the improvement to large collision systems are quite limited. The TBW4 fits are found to work significantly better especially in pp and pPb collisions while χ 2 /nDoF values in central XeXe and PbPb collisions shown in Fig. 3 are quite similar to those obtained in TBW fits presented in Fig. 1. The improvement in small systems mainly comes from a better description to the proton p T distributions as is shown in Fig. 4. This is expected due to the inclusion of the baryon grouped non-extensive parameter q B . This improvement indicates the importance of baryon number in the fragmentation process of small systems. However, it is found that the pull distributions are not flattened to zero especially for protons in large systems even with the independent q B considered. There is a possibility that in central AA collisions, nuclear medium modifications like shadowing effects and parton energy loss effects can be important to understand this discrepancy. From Fig. 3, it is also found that the extracted freeze-out temperature is smaller than that in the default TBW fits for pp and pPb collisions. The non-extensive parameter for meson is usually larger than that for baryon except in the central AA collisions when the off equilibrium effects are vanishing. These two non-extensive parameters q M and q B become similar in central XeXe and PbPb collisions.

B. Kinetic freeze-out parameters
The model parameters extracted from the TBW and TBW4 fits dependent on the event multiplicity are studied in this section. We present the results of kinetic freeze-out temperature T , average radial flow velocity in the transverse plane β and non-extensive parameter q extracted from the TBW analysis varying with the mid-rapidity charged particle multiplicity in Fig. 5. The charged multiplicity is quantified with the average charged particle number within |η| < 0.5 per unit pseudorapidity dN/dη in each event class. The black circles, cyan crosses, red stars, green diamonds, blue triangles and magenta squares represent the results of pp 7TeV, pp 13TeV, pPb 5.02 TeV, XeXe 5.44 TeV, PbPb 2.76 TeV and PbPb 5.02 TeV, respectively. It is shown in Fig. 5(a) that the multiplicity dependent flow velocity can be divided into two groups, one group contains mainly large symmetric systems consisting of XeXe 5.44 TeV, PbPb 2.76 TeV and PbPb 5.02 TeV, the other group includes the small systems ranging from pp 7 TeV, pp 13 TeV to pPb 5.02 TeV. The average flow velocity grows with the multiplicity rapidly in the low dN/dη region but saturates at about 0.6 when the average charge number density dN/dη approaches 1000. The onset of non-zero radial flow starts from dN/dη around 10 to 20 in two groups. The branching of small system flow velocity with respect to the large system flow at the same dN/dη can be understood in the context of parton density difference in the transverse plane. Higher initial densities in pp and pA collisions might produce larger pressure against the surrounding environment [45,46].
The non-extensive parameter is shown in Fig. 5(b) as (q − 1) versus the average charge particle density. It is found that the non-extensive parameter initially increases with the charge particle density and then drops when dN/dη becomes larger than 10. The non-extensive parameters from different collision systems seem to scale together. Deviations arise in the region of dN/dη between 10 and 50. The non-extensive parameter in pPb collisions in the high multiplicity regime follows the AA collision trend, unlike the radial flow velocity which follows the trend of pp collisions. On the other hand, the extracted kinetic freeze-out temperature T shows the opposite dependence on the event multiplicity compared to the non-extensive parameter (q − 1). As is displayed in Fig. 5(c), a slight decreasing within dN/dη < 10 and a mild growth of the kinetic freeze-out temperatures dependent on the charge particle density        in the region of dN/dη > 40 is observed. In both nonextensive parameter and freeze-out temperature comparison, the high multiplicity pp 13 TeV results are found to slightly deviate from the global trend. It is also interesting to see that the radial flow velocity in high multiplicity pPb collisions is similar to pp collisions, while the nonextensive parameter in the same events is following the AA trend.
In the end, we also present the effective temperature of the system defined as Fig. 5(d). The effective temperature incorporating the radial flow velocity seems to be insensitive to the collision system but only relies on the charge multiplicity. A turning point in the T ef f distribution is seen at dN/dη ∼ 10 similar to the non-extensive parameter distribution, above which the temperature increases monotonously up to dN/dη ∼ 10 3 in a unified way across all different systems ranging from pp collisions to AA collisions.
In the framework of non-equilibrium statistics, the temperature and the flow velocity can be related to viscosity with the linear or quadratic dependence on the non-extensive parameter (q − 1). The β vs (q − 1) and T vs (q − 1) distributions from TBW fits are shown in Fig. 6. The non-zero radial flow velocity and (q − 1) can be roughly described by a universal quadratic dependence, with the exception that high multiplicity pp collisions are more likely to have larger non-extensive parameter instead of having vanishing (q − 1) value as seen in central AA collisions. The kinetic freeze-out temperature T vs (q − 1) also shows little sensitivity to the collision system, while a slightly diverged branch consisting of XeXe collisions and pPb collisions is observed. the large symmetric systems is possibly coming from the fluctuation effects in the initial conditions whose initial spatial geometries can vary significantly due to the nonspherical shape of Xe nuclei and the asymmetric colliding nuclei involved in the pPb collisions. The potentially produced hot spots region during the fluctuations may lead longer life time of the evolving medium and a lower freeze out temperature.
The extracted freeze-out parameters are also studied in the case with two independent non-extensive parameters q M for meson and q B for baryon. Compared to the results shown in Fig. 5, much better scaling features are observed in Fig. 7 for TBW4 fits. The radial flow velocity distribution is still divided into two groups with a quite early onset of non-zero radial flow in dN/dη as shown by Fig. 7(a). The non-extensive parameters q B for the baryon are represented by the hollow markers and compared to q M in Fig. 7(b). q M is found to be larger than q B especially in the low multiplicity regime. The separated meson and baryon non-extensive parameters seem to automatically converge in the high multiplicity region. This indicates that the baryon number is mostly important for the hadronization process in the small systems. For large systems, the non-equilibrium effects become negligible and the non-extensive parameter turns out to be independent of the hadron species. The radial flow velocity and the non-extensive parameter in high multiplicity pPb collisions show the same diverged features similar to the finding in Fig. 5, following the pp collision behavior in β and the AA collision behavior in (q − 1). This deviation arises since the transverse overlap area in small systems is much smaller than that in large systems at the same event multiplicities [47]. The initial entropy density, represented by the particle density per unit transverse overlap area, in small systems is thus significantly enhanced, which leads to the stronger collective expansion in high multiplicity pPb collisions in contrast to that in peripheral PbPb collisions. The extracted radial flow velocity is expected to scale with the charged particle density per unit transverse overlap area across the large and small systems in a more coherent way [48].
A two staged temperature dependence varying with dN/dη is found in Fig. 7(c) and Fig. 7(d) similar to the results shown in Fig. 5(c) and Fig. 5(d). The deviation of the temperature from PbPb collisions to XeXe collisions is less important in TBW4 fits than that in the TBW case. Without the initial decreasing part, T ef f in TBW4 grows all the way from low to high multiplicity only with different slopes in the two regions separated by dN/dη ∼ 20.
It is noteworthy that a turning behavior exists in the multiplicity dependence of the freeze-out parameters β , q and T around dN/dη ∼ 10 to 15 shown both in Fig. 5 and Fig. 7. Considering that the hadronic interactions are believed to be dominant in the low multiplicity events and the parton interactions become more important in the high multiplicity limit, the region in between may suffer from the mixture of the two types of contributions. The emergence of this feature suggests that the dominant physics mechanism dictating the evolution of the system changes from the hadron gas rescatterings to the deconfined quark gluon matter interactions. This turning behavior can be regarded as a signature that the quark and gluon degrees of freedom begin to take over in the high multiplicity region.
Considering the TBW4 fits generally give better description to the experimental data, we also examine the β vs (q − 1) and T vs (q − 1) in Fig. 8 from TBW4 analysis. The large system β vs (q − 1) relationship can be roughly described by a quadratic fit, while more complicated structures can be found for the smaller systems as displayed in Fig. 8(a) and Fig. 8(b). Similar to the findings in Fig. 7(b), central pp collisions have large q and sizable flow, while large collision systems tend to have smaller q and large β in central events. The parameters of peripheral pPb events are close to the pp case but become similar to the AA case in high multiplicity pPb collisions. T vs (q − 1) relationships shown in Fig. 8(c) and Fig. 8(d) are found to be universal for all collision systems. Central pp collisions approach large q and small T while central AA collisions are found to have large T and small q. It is interesting to see that both the AA system and the pp system have similar model parameter values in the peripheral collisions. These two types of collision systems are following two distinguished evolution curve on the β vs (q − 1) plane and moving towards opposite direction along a universal curve in the T vs (q − 1) space. The difference between the large systems and the small systems in the β vs (q − 1) space can be partly understood after considering the striking β gap between the two types of systems shown in Fig. 5(a) and Fig. 7(a). Additionally, the initial spatial fluctuation effects can be also important to account for this behavior. The initial geometries of large systems are supposed to be controlled by the average shape of the overlap region while the sub-nucleon level fluctuation effects arise in the determination of initial densities for the high multiplicity small systems [49,50]. The sizable non-equilibrium effects induced by the strong fluctuations in central pp collisions lead to increasing non-extensive parameter, unlike the central AA collisions in which q is vanishing as the system approaches equilibrium. As the dynamical fluctuations in freeze-out temperature induced by the system size effects are encoded in the non-extensive parameter [12,51,52], the variations in temperature are compensated by the corresponding non-extensive parameters, universal scalings are observed in the T vs (q − 1) space. It is expected that the pp collisions and pPb collisions with extremely high multiplicities will approach the environment created in AA collisions, indicated by the behavior of β and (q − 1) in Fig. 5 and Fig. 7. Examining whether the turning behavior shown in Fig. 8(a) and (b) exists in the ultra-central pp and pPb collisions is an interesting analysis with the future experiment data.

IV. SUMMARY
In this paper, we have used the Tsallis Blast-Wave model with and without independent baryon nonextensive parameters to fit the transverse momentum spectra of charged pion, kaon and protons produced at mid-rapidity in pp collisions at √ s = 7 and 13 TeV, pPb collisions at √ s NN = 5.02 TeV, XeXe collisions at √ s NN = 5.44 TeV, and in PbPb collisions at √ s NN = 2.76 TeV and 5.02 TeV at the LHC to extract kinetic freeze-out parameters. It is found that the introduction of q B for the baryon in the TBW4 fits improves the description of identified hadron spectra especially in pp, pA and peripheral AA collisions. The multiplicity dependence of the freeze-out properties are examined across the the largely varied collision systems. A general universal scaling of the freeze-out parameters can be observed especially in the Tsallis analysis with separate baryon non-extensive parameter. The radial flow velocities obtained from both the TBW fits and the TBW4 fits are divided into two categories consisting of small systems and large symmetric systems, respectively. The temperatures from different collision systems are scaling with the charge multiplicity in a unified way. Similar universality is also found in the nonextensive parameter distributions with some minor deviations in the most central pp collisions. By investigating the correlation of radial flow velocity β vs (q − 1) and kinetic freeze-out temperature T vs (q − 1) in the TBW4 fits with independent q B assumptions, the peripheral AA collisions and pp collisions are found to be quite similar while the two types of systems are moving towards completely different direction going to central collisions in the parameter space. The asymmetric pPb collisions are close to the pp collisions in the low multiplicity region but become AA like in the high multiplicity region. A transitional behavior can be found in the evolution of the pPb collision system.
The universality of the freeze-out properties revealed in this work suggests that the Tsallis Blast-Wave analysis can be applicable to various collision systems covering a wide range of event multiplicities. This feature also indicates that the evolution properties of the small size collision systems with high multiplicities and the large systems at the LHC energy can be driven by a unified physics mechanism. A parton evolution stage with sizable collective motion might exist even in the small collision systems. Extending the study to the pp collision events with even higher charge particle density [53] or to the intermediate size collision systems like Oxygen-Oxygen collisions [54] in the future experiments can be important to further understand the universality of freeze-out properties observed in different collision systems at the LHC energies.