Analysis of Midrapidity p T Distributions of Identiﬁed Charged Particles in Pb + Pb Collisions at √ s nn = 5.02 TeV Using Tsallis Distribution with Embedded Transverse Flow

: The midrapidity transverse momentum distributions of the charged pions, kaons, protons, and antiprotons, measured by ALICE Collaboration at ten centrality classes of Pb + Pb collisions at √ s nn = 5.02 TeV in the Large Hadron Collider (LHC, CERN, Switzerland), are successfully analyzed using combined minimum χ 2 fits with a thermodynamically non-consistent, as well as thermodynamically consistent, Tsallis function with transverse flow. The extracted non-extensivity parameter q decreases systematically for all considered particle species with increasing Pb + Pb collision centrality, suggesting an increase in the degree of system thermalization with an increase in collision centrality. The results for q suggest quite a large degree of thermalization of quark–gluon plasma (QGP) created in central Pb + Pb collisions at √ s nn = 5.02 TeV with the average number of participant nucleons (cid:10) N part (cid:11) > 160. The obtained significantly different growth rates of transverse flow velocity, (cid:104) β T (cid:105) , in regions (cid:10) N part (cid:11) < 71 ± 7 and (cid:10) N part (cid:11) > 71 ± 7 with the temperature parameter T 0 remaining constant within uncertainties in region (cid:10) N part (cid:11) > 71 ± 7 probably indicates that (cid:10) N part (cid:11) ≈ 71 ± 7 (corresponding to (cid:104) dN ch / d η (cid:105) ≈ 251 ± 20) is a threshold border value for a crossover transition from a dense hadronic state to the QGP phase (or mixed phase of QGP and hadrons) in Pb + Pb collisions at √ s nn = 5.02 TeV. The threshold border value for transverse flow velocity (cid:104) β T (cid:105) ≈ 0.46 ± 0.03 (corresponding to (cid:10) N part (cid:11) ≈ 71 ± 7), estimated by us in Pb + Pb collisions at √ s nn = 5.02 TeV, agrees well with the corresponding border value (cid:104) β T (cid:105) ≈ 0.44 ± 0.02, recently obtained in Xe + Xe collisions at √ s nn = 5.44 TeV, and with almost constant (cid:104) β T (cid:105) values extracted earlier in the Beam Energy Scan (BES) program of the Relativistic Heavy-Ion Collider (RHIC, Brookhaven, GA, USA) in central Au + Au collisions in the √ s nn = 7.7 − 39 GeV energy range, where the threshold for QGP production is achieved. The correlations between extracted T 0 and (cid:104) β T (cid:105) parameters are found to be greatly different in regions (cid:104) β T (cid:105) < 0.46 and (cid:104) β T (cid:105) > 0.46, which further supports our result obtained for the threshold border value in Pb + Pb collisions at √ s nn = 5.02 TeV.


Introduction
It is one of the main goals of modern experiments conducted using the RHIC (Relativistic Heavy-Ion Collider, Brookhaven, NY, USA) and LHC (Large Hadron Collider, CERN, Meyrin, Switzerland) to produce in collisions of high-energy heavy ions and investigate in detail quark-gluon plasma (QGP), the plasma state of almost-free quarks and gluons. Initially produced QGP expands hydrodynamically at an extremely high speed, reaching, one after the other, the chemical and kinetic freeze-out stages, respectively, which fixes the hadron abundances (yields) and final momenta of hadrons. Because of the extremely short lifespan of QGP, it is impossible to detect and measure its properties directly. Therefore, to extract important information on QGP and its evolution with changing collision energy and centrality, one has to analyze the yields, properties, and transverse momentum or/and (pseudo)rapidity spectra of the final particles with the help of efficient theoretical and phenomenological models. Pions, kaons, and (anti)protons account for the predominantly largest part of the final particles produced in high-energy collisions in the RHIC and LHC. Therefore, the production of these particles, consisting of light (u, d, and s) quarks, has been studied quite extensively  to deduce valuable information on the evolution and properties of produced hot and dense matter, and mechanisms of particle production at high-energy collisions.
Different modifications of the so-called Tsallis distribution function, based on nonextensive Tsallis statistics, have become efficient and standard functions for describing the transverse momentum distributions of particles in proton + proton collisions at high energies [71][72][73][74][75][76][77][78][79][80][81][82][83]. Compared to other power-law distributions, the Tsallis distribution function has a clear advantage expressed by its direct connection to thermodynamics through entropy [80]. In addition to the effective temperature parameter, the non-extensivity index q of the Tsallis function represents an important parameter, which shows the degree of deviation of particle transverse momentum distribution from the exponential Boltzmann-Gibbs distribution. Moreover, parameter q is said to measure the degree of non-equilibrium or that of the non-thermalization of a system [84]. The importance of the q parameter and its deep physical meaning, with a direct connection to thermodynamics, and the physical significance of q in non-extensive statistical mechanics have been proven once again in the recent work [71] of Tsallis. In Refs. [23,49,[62][63][64]85], the transverse flow velocity was embedded using a simple Lorentz transformation into a QCD-inspired (power law) Hagedorn function to successfully reproduce the measured longer ranges of p T spectra of the final particles in heavy-ion as well as proton + proton collisions at high energies. In a recent study [24], using the simple Lorentz transformation, we incorporated for the first time the transverse expansion velocity into both the simple, thermodynamically non-consistent, and thermodynamically consistent Tsallis functions to successfully analyze the centrality dependencies of p T distributions of the charged pions, kaons, and (anti)protons in Xe + Xe collisions at √ s nn = 5.44 TeV, measured by ALICE Collaboration in Ref. [12]. Compared to other hydrodynamics-based blast-wave models, the non-consistent Tsallis function with a transverse flow, as well as the thermodynamically consistent Tsallis function with a transverse flow, obtained for the first time in Ref. [24] and used in the present analysis, have fewer parameters, and each of them carries a physical meaning [24]. As indicated in Refs. [1,49,57,63,64], it is impossible to assign any physical meaning to collectivity parameters, such as transverse flow velocity and kinetic freeze-out temperature, obtained from separate model fits to p T distribution of single-particle species. On the other hand, the simultaneous combined model fits to p T distributions of all analyzed particle species in a given collision system can produce the physically meaningful collectivity parameters of the analyzed system [1,20,23,24,49,57,[62][63][64]85]. The combined (global) model fits permit us to compare various collision systems using a few parameters [1,20,23,24,49,57,[62][63][64]85]. The present study aims to investigate the midrapidity (mid-y) p T distributions of identified charged particles, measured by the ALICE Collaboration [57] at different centralities of Pb + Pb collisions at √ s nn = 5.02 TeV, applying combined model fits with both the non-consistent as well as thermodynamically consistent Tsallis functions with transverse flow, introduced and used for the first time in Ref. [24] to successfully study the collision centrality dependencies of mid-y p T distributions of the identified charged particles in Xe + Xe collisions at √ s nn = 5.44 TeV. The paper is organized as follows: the information about the experimental data and theoretical models is presented in Section 2; the analysis and results are presented in Section 3; and Section 4 summarizes the results and conclusions of the present analysis.

Experimental Data and Models
In the present study, we analyzed the experimental data obtained by ALICE Collaboration [57] at the LHC using the central barrel of the ALICE detector having full azimuthal coverage in the mid-(pseudo)rapidity |η| < 0.8 region [86]. In order to obtain transverse momentum distributions of the primary charged particles, the raw p T distributions were corrected by ALICE Collaboration for misidentification probability, acceptance, reconstruction, and tracking efficiencies, and contaminations mostly due to weak decays of Λ and Σ + (affecting (anti)proton spectra), and of K s 0 (affecting spectra of charged pions), and interactions with the material were computed and subtracted from the raw spectra. The p T spectra were extracted for 0-5%, 5-10%, 10-20%, 20-30%, 30-40%, 40-50%, 50-60%, 60-70%, 70-80%, and 80-90% collision centralities [57]. The average values of the number of participant nucleons ( N part ) and charged-particle (pseudo-rapidity) multiplicity density ( dN ch /dη ) were obtained [57,87] by ALICE Collaboration using Glauber-Monte Carlo model calculations for the analyzed centrality groups of Pb + Pb collisions at √ s nn = 5.02 TeV, which are presented in Table 1. Table 1. The average number of participant nucleons and mean charged-particle multiplicity densities [57,87] in the analyzed groups of centralities of Pb + Pb collisions at √ s nn = 5.02 TeV. There are various forms [62,75,76,[80][81][82][83]88,89] of the Tsallis distribution function. They all demonstrate quite good quality fits of the experimental p T spectra of particles originating from high-energy collisions. One of the simplest versions [82] of the Tsallis distribution is presented at mid-y (<y> = 0) by the following expression:

Centr
where C is the fitting constant, N ev is the total number of inelastic collision events, m T = p 2 T + m 2 0 is the transverse energy, m 0 is the rest mass of a hadron, T represents the effective temperature, and q is the so-called non-extensivity parameter. This function in Equation (1) is called the simple, thermodynamically non-consistent Tsallis distribution in the present study. Parameter q is important, and it shows the deviation of transverse momentum distributions from the exponential Boltzmann-Gibbs distribution. This parameter is also said to measure the non-equilibrium level [84]. At q = 1, the Tsallis distribution reduces to the exponential (equilibrated) Boltzmann-Gibbs distribution.
The following form [75,[80][81][82] of the Tsallis function at mid-y (<y> = 0) obtained at zero chemical potential µ (µ approaches zero and is negligible at high-collision energies at the LHC) results in consistent thermodynamics for the pressure, temperature, particle number, energy, and entropy densities: which is referred to as a consistent Tsallis function or Tsallis function with thermodynamical consistence in the present analysis. The thermodynamical consistence of this function has been proven in Ref. [75]. Here, C q is the fitting constant, which is related linearly to the system's volume (V) [75,80] through the expression C q = gV/(2π) 3 , where g is the degeneracy factor equal to 2, 3, 4 for protons, pions, and kaons, respectively.
The T values extracted from fits by Tsallis functions in Equations (1) and (2) denote the effective temperature, which contains the effects of both the thermal motion and collective flow (expansion). The transverse expansion velocity should be embedded into the Tsallis distribution function in order to disentangle these two effects.
In Refs. [23,49,[62][63][64]85], the transverse flow was incorporated into a QCD-inspired Hagedorn function, 2πN ev p T dp T dy = C n 1 + m T p 0 −n , by setting p 0 = nT 0 and applying the simple Lorentz transformation m T → γ T ·(m T − p T β T ) to create the final expression: where C n is the fitting constant, γ T = 1/ 1 − β T 2 , < β T > is the average transverse flow velocity, T 0 is the parameter estimating the kinetic freeze-out temperature, and n is the free parameter. The function in Equation (3) is called the Hagedorn function with an (embedded) transverse flow. The important substitution, m T → γ T ·(m T − p T β T ) , used to derive Equation (3) represents the Lorentz transformation into a system co-moving with the average (common) transverse flow velocity, < β T >, of particles if such a flow of particles exists in the co-moving frame [62]. The Hagedorn function with transverse flow, presented in Equation (3), can describe quite well the p T spectra of particles produced in proton + proton and heavy-ion collisions in the RHIC and LHC in Refs. [23,49,[62][63][64]85], allowing for the clear physical interpretation of the extracted parameters.
In comparison to Refs. [23,49,[62][63][64]85], where the transverse expansion velocity was incorporated into the Hagedorn function resulting in the final Equation (3), in the present study, we embedded the transverse flow velocity into both the simple Tsallis function (Equation (1)) and Tsallis function with thermodynamical consistence (Equation (2)), as was performed for the first time in Ref. [24]. To describe mid-y p T distributions, d 2 N N ev dp T dy , of the identified charged particles in different centrality classes of Pb + Pb collisions at √ s nn = 5.02 TeV, we embedded the transverse expansion velocity into the simple Tsallis function in Equation (1) by replacing m T → γ T ·(m T − p T β T ) [24]: which is called the simple or non-consistent Tsallis function with a transverse flow in the present study. One can observe that Equations (3) and (4) have similar mathematical structures and are almost equivalent, with the only difference being in parameters n and q. Similarly, we inserted the transverse flow into the Tsallis function with thermodynamical consistence in Equation (2) by performing the m T → γ T ·(m T − p T β T ) transformation to obtain [24] which is called a thermodynamically consistent Tsallis function with a transverse flow in the present study. It is necessary to note that the functions in Equations (4) and (5) were derived for the first time and used successfully in Ref. [24] to describe the mid-y p T distributions of the identified charged particles in different centrality classes of Xe + Xe collisions at √ s nn = 5.44 TeV.
In the present study, we conducted simultaneous (combined) fits by the model functions, presented in Equations (4) and (5), of transverse momentum distributions of the studied particle species in each of the ten centrality groups employing the nonlinear curve fitting of Origin 9.1 Graphing and Data Analysis Software. T 0 and β T were kept as the common (shared) parameters for all studied particle species during the combined fitting procedures.
The experimental data points in the figures of the present study show the combined (added) statistical and systematic errors. The combined errors are determined predominantly by the systematic uncertainties, which are much greater than the statistical ones. More details on the calculation of the systematic errors in the measured p T distributions can be found in Ref. [57]. The minimum χ 2 fitting procedures were made taking into account the combined statistical and systematic errors as the weights (1/(error) 2 ) for the data points. The region p T < 0.5 GeV/c in pion p T distributions, which has a significant contribution from decays of baryon resonances, was excluded from fitting procedures, similar to Refs. [1,23,24,49,57,63,64,85]. In the present study, we used the same p T intervals for combined fits with the model functions in Equations (4) and (5) as in Ref. [24]: [0.5-5.0] GeV/c for π + + π − , [0. [2][3][4][5].0] GeV/c for K + + K − , and [0. [3][4][5].0] GeV/c for p + p. By using the identical p T intervals and model functions, and the same analyzed particle species as those in Ref. [24], we could directly compare the results obtained in the present study for Pb + Pb collisions at √ s nn = 5.02 TeV with those extracted recently in Ref. [24] for Xe + Xe collisions at √ s nn = 5.44 TeV.

Analysis and Results
The results of combined minimum χ 2 fits using the thermodynamically consistent Tsallis function with a transverse flow (Equation (5)) of midrapidity transverse momentum distributions of the charged pions, kaons, protons, and antiprotons in ten groups of centrality of Pb + Pb collisions at √ s nn = 5.02 TeV are presented in Table 2. The corresponding combined minimum χ 2 fit curves with a consistent Tsallis function with a transverse flow in four different centrality classes of Pb + Pb collisions are presented in Figure 1. Table 2. The results of combined minimum χ 2 fits with thermodynamically consistent Tsallis function with transverse flow (Equation (5)) of p T spectra of particles in Pb + Pb collisions at √ s nn = 5.02 TeV.
n.d.f. denotes the number of degrees of freedom.  Table 3 presents the parameters obtained from combined minimum χ 2 fits using the thermodynamically non-consistent Tsallis function with a transverse flow (Equation (4)) of mid-y p T distributions of the studied particle species in ten groups of centrality of Pb + Pb collisions at √ s nn = 5.02 TeV. The corresponding combined minimum χ 2 fit curves by the non-consistent Tsallis function with transverse flow in four different centrality classes of Pb + Pb collisions are illustrated in Figure 2. Table 3. The results of combined minimum χ 2 fits with thermodynamically non-consistent Tsallis function with transverse flow (Equation (4)) of p T spectra of particles in Pb + Pb collisions at √ s nn = 5.02 TeV.     Tables 2 and 3 show that the combined fits with both consistent and non-consistent Tsallis functions with transverse flows reproduce quite well the midrapidity transverse momentum distributions of the studied particle species in ten centrality classes of Pb + Pb collisions at √ s nn = 5.02 TeV. It can be observed in Tables 2 and 3 that the fits by the simple non-consistent Tsallis function with a transverse flow result in slightly lower values of χ 2 /n.d.f. than those by the consistent Tsallis function with transverse flow. Though the quality of the fits by the function in Equation (4) proved to be slightly better than those by the function in Equation (5) from the mathematical point of view, from a physics perspective it is still advantageous to use the function in Equation (5), based on the Tsallis distribution with thermodynamical consistence. This is because the consistent Tsallis distribution satisfies all the thermodynamic relations [75,80,81] (those for the pressure, temperature, entropy, and particle densities, etc.) resulting from the equations of the first and second laws of thermodynamics. The clear advantage is, therefore, that temperature parameter T in the Tsallis distribution with thermodynamical consistence [75] can be considered as true temperature in the thermodynamic sense.  Tables 2 and 3 show that the combined fits with both consistent and non-consistent Tsallis functions with transverse flows reproduce quite well the midrapidity transverse momentum distributions of the studied particle species in ten centrality classes of Pb + Pb collisions at = 5.02 TeV. It can be observed in Tables 2 and 3 that the fits by the simple non-consistent Tsallis function with a transverse flow result in slightly lower values of /n.d.f. than those by the consistent Tsallis function with transverse flow. Though the quality of the fits by the function in Equation (4) proved to be slightly better than those by the function in Equation (5) from the mathematical point    Tables 2 and 3. In addition, the obtained β T versus dN ch /dη dependence is also presented. The analysis of β T and T 0 behaviors and their dependencies on collision centrality and multiplicity of particles is important [23,24,49,66] because these parameters are also related to mapping the QCD phase diagram; however, the chemical freeze-out temperature is commonly used in such phase diagrams.   (5)), presented in Table  2, for the charged pions (•), kaons (▲), and protons and antiprotons (▪). The results (Table 3)   ). The results (Table 3) (8)) and two-power (Equations (7) and (9)) functions, respectively, of the β T versus N part and β T versus dN ch /dη dependencies, respectively, extracted using consistent Tsallis function with transverse flow. For better visibility, the data (open symbols), extracted using non-consistent Tsallis function with transverse flow (Equation (4)), are slightly shifted along the positive directions of N part and dN ch /dη axes. Figure 3a,b demonstrate that parameter T 0 as well as β T , obtained using both consistent and non-consistent Tsallis functions with transverse flows, similarly depend on N part . Figure 3a shows that T 0 values obtained from combined fits by the consistent Tsallis function with a transverse flow (Equation (5)) are systematically smaller as compared to those from fits by a non-consistent Tsallis function with transverse flow (Equation (4)). This can be explained by the extra γ T ·(m T − p T β T ) factor in Equation (5), which is absent in Equation (4).
As can be observed in Figure 3a,b, the N part dependencies of parameters T 0 and β T in non-central Pb + Pb collisions in region N part < 71 ± 7 differ substantially from those in semi-central and central Pb + Pb collisions in region N part > 71 ± 7. The border value N part ≈ 71 ± 7 is estimated as the middle between the 4th and 5th points on the N part axis (corresponding to 50-60% and 40-50% centrality in Table 1) in Figure 3a,b. Figure 3a shows that parameter T 0 on the whole decreases significantly in region N part < 71 ± 7, then T 0 becomes quite stable being constant within fit uncertainties in region N part > 71 ± 7. As observed from Figure 3b As observed from Figure 3b, the parameter β T demonstrates a fairly high growth rate at N part < 71 ± 7 and substantially lower rate of increase at N part > 71 ± 7. To prove this observation, we fitted β T versus N part dependence in Figure 3b, obtained using the consistent Tsallis function with a transverse flow, with a single-power function as well as two-power functions: where u(t) is the Heaviside function (u(t) = 0 if t < 0, and u(t) = 1 if t > 0), A 1 and A 2 are fitting (normalization) constants, and α 1 and α 2 are exponent parameters. As can be observed in Figure 3c, the parameter β T demonstrates quite a high growth rate at dN ch /dη < 251 ± 20 and a substantially smaller rate of increase at dN ch /dη > 251 ± 20. The threshold border value dN ch /dη ≈ 251 ± 20, corresponding to border value N part ≈ 71 ± 7, was estimated as the middle point between dN ch /dη values for 50-60% and 40-50% collision centralities in Table 1. Similarly to Figure 3b, we fitted the β T versus dN ch /dη dependence in Figure 3c, obtained using consistent the Tsallis function with transverse flow, with a single-power function β T = A· dN ch /dη α ·(A-fitting constant, α-exponent parameter) as well as a two-power function: Though, for the convenience, we used the same symbols for the parameters in Equations (6) and (8) (and in Equations (7) and (9)), they should be considered as different parameters in equations with different arguments of N part and dN ch /dη . The parameters obtained from minimum χ 2 fits of β T versus N part ( β T versus dN ch /dη ) dependence in Figure 3b (Figure 3c) by the functions in Equations (6) and (7) (Equations (8) and (9)) are presented in Table 4 (Table 5). As can be observed in Figure 3b ( Figure 3c) and χ 2 /n.d.f. values in Table 4 (Table 5), the single-power function cannot satisfactorily describe the β T versus N part ( β T versus dN ch /dη ) dependence. At the same time, the two-power function with two different exponent parameters in regions N part < 71 ± 7 ( dN ch /dη < 251 ± 20) and N part > 71 ± 7 ( dN ch /dη > 251 ± 20) reproduces these dependencies in Figure 3b,c very well. Significantly differing growth rates of β T in regions N part < 71 ± 7 ( dN ch /dη < 251 ± 20) and N part > 71 ± 7 ( dN ch /dη > 251 ± 20) with T 0 staying constant within uncertainties in region N part > 71 ± 7 ( dN ch /dη > 251 ± 20) possibly suggests that N part ≈ 71 ± 7 ( dN ch /dη ≈ 251 ± 20) is a threshold border value for a crossover transition from a dense hadronic (nucleon) state to the QGP phase (or mixed phase of QGP and hadrons) in Pb + Pb collisions at √ s nn = 5.02 TeV. We can estimate the average transverse flow velocity, β T , corresponding to a possible border value N part ≈ 71 ± 7 ( dN ch /dη ≈ 251 ± 20) for a crossover transition to the QGP state (or mixed phase of QGP and hadrons) in Pb + Pb collisions. The threshold border value β T ≈ 0.46 ± 0.03, corresponding to N part ≈ 71 ± 7 ( dN ch /dη ≈ 251 ± 20), was estimated as the middle between β T values for 50-60% and 40-50% collision centralities in Table 2, obtained from combined fits with the function in Equation (5) of transverse momentum distributions of the considered particle species in Pb + Pb collisions at √ s nn = 5.02 TeV. Table 4. The results of minimum χ 2 fits with the functions in Equations (6) and (7) of β T versus N part dependence in Figure 3b, obtained using consistent Tsallis function with transverse flow (Equation (5)).

Equation (6)
A = 0.143 ± 0.025 α = 0.252 ± 0.033 159/8 Equation (7) A 1 = 0.052 ± 0.007 α 1 = 0.532 ± 0.034 A 2 = 0.282 ± 0.022 α 2 = 0.129 ± 0.014 6/6 Table 5. The results of minimum χ 2 fits with the functions in Equations (8) and (9) of β T versus dN ch /dη dependence in Figure 3c, obtained using consistent Tsallis function with transverse flow (Equation (5)). It is important to mention that similar results were obtained in our recent work [24] from the analysis of transverse momentum distributions of identified charged particles in different centrality classes of Xe + Xe collisions at √ s nn = 5.44 TeV. The results, similar to those presented in Figure 3, but obtained in Ref. [24] for N part dependencies of the extracted parameters T 0 , β T , and q in Xe + Xe collisions at √ s nn = 5.44 TeV, are demonstrated in Figure 4. The average transverse flow velocity, β T , presented quite a large growth rate in region N part < 44 ± 5 and a significantly smaller growth rate in range N part > 44 ± 5 with increasing N part in Xe + Xe collisions at √ s nn = 5.44 TeV [24].

Fitting Function
Interestingly, the probable crossover transition border to the QGP state (or mixed phase of QGP and hadrons) was estimated to be between 50-60% and 40-50% collision centralities coinciding for both Xe + Xe collisions at √ s nn = 5.44 TeV in Ref. [24] and Pb + Pb collisions at √ s nn = 5.02 TeV in the present analysis. The significantly different growth rates of β T in regions N part < 44 ± 5 and N part > 44 ± 5 with parameter T 0 stabilizing and becoming practically constant in range N part > 44 ± 5 were observed in Xe + Xe collisions [24], indicating that N part ≈ 44 ± 5 (corresponding to dN ch /dη ≈ 158 ± 20) could be a threshold border value for a crossover transition from a dense hadronic state to the QGP phase (or mixed phase of QGP and hadrons) in Xe + Xe collisions at √ s nn = 5.44 TeV.
It is interesting to note that the threshold border value β T ≈ 0.44 ± 0.02 (corresponding to N part ≈ 44 ± 5 and dN ch /dη ≈ 158 ± 20), recently estimated in Xe + Xe collisions at √ s nn = 5.44 TeV in Ref. [24], agreed within uncertainties with the corresponding fi T ≈ 0.46 ± 0.03 obtained for Pb + Pb collisions at √ s nn = 5.02 TeV in the present analysis, and with nearly constant β T values extracted from the BES program at the RHIC [20] in central Au + Au collisions in the √ s nn = 7.7−39 GeV energy range, where the threshold for QGP production is surely reached [20,[90][91][92][93][94].
It is worth presenting in more detail the results of Ref. [20], where the chemical and kinetic freeze-out temperatures and average transverse flow velocities extracted in central (0-5%) Au + Au collisions at BES energies from √ s nn = 7.7 to 39 GeV are compared with those obtained in central (0-5% and 0-7% centrality) heavy-ion collisions at AGS, SPS, top RHIC, and up to LHC energies (see Figure 38 in Ref. [20]). The extracted β T increased at low √ s nn < 7.7 GeV energies and remained nearly constant across BES energies in central Au + Au collisions, changing from 0.44 ± 0.04 at √ s nn = 7.7 GeV to 0.49 ± 0.04 at √ s nn = 39 GeV [20]. After it, β T gradually and steadily increased to LHC energies [20]. The obtained chemical freeze-out temperature presented an increase from √ s nn = 7.7 to 19.6 GeV, then remained constant within the errors across higher BES, top RHIC, and up to LHC energies [20]. The chemical freeze-out temperature proved to be ≈ 160 ± 5 MeV in central heavy-ion collisions in the √ s nn ≥ 19.6 GeV range [20], agreeing with the critical deconfinement phase transition temperature, T c ≈ 150-170 MeV, coming from lattice QCD calculations [36,37]. The yield of pions is characterized by a linear increase with increasing √ s nn at lower collision energies, and it shows a kink structure at √ s nn ≈ 19.6 GeV in central Au + Au collisions at the BES [20]. This has been interpreted as a clear change in the mechanism of particle production at √ s nn ≈ 19.6 GeV [20].
As observed in Figure 3d, the parameter of non-extensivity q for the charged pions, kaons, and (anti)protons, resulting from both consistent and non-consistent Tsallis functions with transverse flow, shows similar dependencies on N part . The q values demonstrate a systematic decrease with increasing N part and collision centrality for all studied particle species. A similar result, as observed in Figure 4c, is obtained for q versus N part dependencies for the analyzed particle species in Xe+Xe collisions at √ s nn = 5.44 TeV in Ref. [24]. This observation is in agreement with the fact that the level of thermalization of the system increases with increasing centrality of heavy-ion collisions at high energies. The relation q(pions) ≈ q(kaons) > q(anti(protons)) was satisfied, as observed in Figure 3c, which agrees with the q(mesons) > q(baryons) relation found previously in high-energy collisions in Refs. [23,49,63,82,83,85,95]. On the whole, the gap between q(mesons) and q(baryons) decreases with an increase in Pb + Pb collision centrality, as observed in Figure 3c, and q(mesons) practically coincides within uncertainties with q(baryons) in central collisions with N part > 160. This could indicate quite a large degree of thermalization of QGP produced in central Pb + Pb collisions at √ s nn = 5.02 TeV with N part > 160. It is worth presenting in more detail the results of Ref. [20], where the chemical and kinetic freeze-out temperatures and average transverse flow velocities extracted in central (0-5%) Au + Au collisions at BES energies from = 7.7 to 39 GeV are compared with those obtained in central (0-5% and 0-7% centrality) heavy-ion collisions at AGS, SPS, top RHIC, and up to LHC energies (see Figure 38 in Ref. [20]). The extracted 〈 〉 increased at low < 7.7 GeV energies and remained nearly constant across BES energies in central Au + Au collisions, changing from 0.44 ± 0.04 at = 7.7 GeV to 0.49 ± 0.04 at = 39 GeV [20]. After it, 〈 〉 gradually and steadily increased to LHC energies [20]. The obtained chemical freeze-out temperature presented an increase from = 7.7 to 19.6 GeV, then remained constant within the errors across higher BES, top RHIC, and up to LHC energies [20]. The chemical freeze-out temperature proved to be ≈ 160 ± 5 MeV in central  Figure 3, but obtained in Ref. [24] for N part dependencies of the extracted parameters T 0 (a), β T (b), and q (c) in Xe+Xe collisions at √ s nn = 5.44 TeV. The dashed and solid curves in panel (b) are minimum χ 2 fits by the simple (single)-power and two-power functions, respectively, of the β T versus N part dependence, obtained using thermodynamically consistent Tsallis function with transverse flow (Equation (5)). For better visibility, the data (open symbols) extracted using thermodynamically non-consistent Tsallis function with transverse flow (Equation (4)) are shifted by +3 units along the positive direction of the N part axis.
Quite a significant correlation exists [1,24] between the temperature and flow velocity parameters. To verify it quantitatively, we studied the correlations between the extracted parameters T 0 and β T , presented in Tables 2 and 3. The dependencies of T 0 versus β T parameters extracted from fits with consistent (Equation (5)) and non-consistent (Equation (4)) Tsallis functions with transverse flow and presented in Tables 2 and 3 are illustrated in Figure 5a,b, respectively. The 1-sigma confidence ellipse of the covariance of the parameters T 0 and β T , corresponding to a 68% confidence interval, and the calculated Pearson correlation coefficient, r xy , between parameters T 0 and β T are also presented in the figures. It is necessary to mention that the Pearson correlation coefficient indicates the degree of a linear correlation between two sets of data (two parameter sets), and it can adopt values between −1 and +1. The orientations and shapes of the confidence ellipses and values of r xy in Figure 5a,b show the high degree of anticorrelation (negative correlation) between parameters T 0 and β T for fits by both consistent and non-consistent Tsallis functions with transverse flows. As can be observed in Figure 5a,b, the r xy values are close to −1 in both cases, and they also indicate a noticeably greater degree of negative correlation between parameters T 0 and β T in Table 3 than that in Table 2. tive correlation in region 〈 〉 < 0.46, the correlation between parameters T0 and 〈 〉 is significantly positive in region 〈 〉 > 0.46 in Figures 6b,d. Hence, it can be observed that not only does parameter 〈 〉 show significantly differing growth rates in regions 〈 〉 < 71 ± 7 (〈 / 〉 < 251 ± 20) and 〈 〉 > 71 ± 7 (〈 / 〉 > 251 ± 20), but also the character of correlation between parameters T0 and 〈 〉 differs considerably in the corresponding intervals 〈 〉 < 0.46 and 〈 〉 > 0.46. This further supports our conjecture that 〈 〉 ≈ 71 ± 7 (〈 / 〉 ≈ 251 ± 20) is probably a threshold border value for a crossover transition from a dense hadronic (nucleon) state to the QGP phase (or mixed phase of QGP and hadrons) in Pb + Pb collisions at = 5.02 TeV.  (5)) and presented in Table 2. (b)-Dependence of T0 versus 〈 〉 parameters (•) obtained from fits with non-consistent Tsallis functions with transverse flows (Equation (4)) and presented in Table 3. The 1-sigma confidence ellipse (corresponding to a 68% confidence interval) of the covariance of parameters T0 and 〈 〉 and the calculated Pearson correlation coefficient, rxy, between T0 and 〈 〉 are also shown in the figures.  (5)) and presented in Table 2. (b)-Dependence of T 0 versus β T parameters (•) obtained from fits with non-consistent Tsallis functions with transverse flows (Equation (4)) and presented in Table 3. The 1-sigma confidence ellipse (corresponding to a 68% confidence interval) of the covariance of parameters T 0 and β T and the calculated Pearson correlation coefficient, r xy , between T 0 and β T are also shown in the figures.
We presented above that the extracted parameter β T demonstrates significantly differing growth rates in regions N part < 71±7 and N part > 71 ± 7, and T 0 remains practically constant in region N part > 71 ± 7, which probably indicates that N part ≈ 71 ± 7 ( dN ch /dη ≈ 251±20) is a border value for a crossover transition from a dense hadronic state to the QGP phase (or mixed phase of QGP and hadrons) in Pb + Pb collisions at √ s nn = 5.02 TeV. The border value β T ≈ 0.46 ± 0.03, corresponding to N part ≈ 71±7, was estimated above from the β T data presented in Table 2. As can be observed in Figure 5a, the characters of T 0 versus β T dependence are remarkably different in regions N part < 71 ± 7 and N part > 71 ± 7. The analogous situation can be observed in Figure 5b. To check and quantify this observation, we studied the correlations between parameters T 0 and β T , presented in Figure 5a,b, separately in regions β T < 0.46 and β T > 0. 46. The results of the separate analysis of the data presented in Figure 5a,b (Tables 2 and 3) in region β T < 0.46 are shown in Figure 6a,c, respectively. The corresponding results for the correlation analysis of the same data in region β T > 0.46 are demonstrated in Figure 6b,d, respectively. As can be observed in Figure 6a,c, the negative correlation between T 0 and β T is quite strong in region β T < 0.46. In contrast to the strong negative correlation in region β T < 0.46, the correlation between parameters T 0 and β T is significantly positive in region β T > 0.46 in Figure 6b,d. Hence, it can be observed that not only does parameter β T show significantly differing growth rates in regions N part < 71 ± 7 ( dN ch dη < 251 ± 20) and N part > 71 ± 7 ( dN ch dη > 251 ± 20), but also the character of correlation between parameters T 0 and β T differs considerably in the corresponding intervals β T < 0.46 and β T > 0.46. This further supports our conjecture that N part ≈ 71 ± 7 ( dN ch /dη ≈ 251 ± 20) is probably a threshold border value for a crossover transition from a dense hadronic (nucleon) state to the QGP phase (or mixed phase of QGP and hadrons) in Pb + Pb collisions at √ s nn = 5.02 TeV.  (5)) and presented in Table 2. (b)-Dependence of T0 versus 〈 〉 parameters (•) obtained from fits with non-consistent Tsallis functions with transverse flows (Equation (4)) and presented in Table 3. The 1-sigma confidence ellipse (corresponding to a 68% confidence interval) of the covariance of parameters T0 and 〈 〉 and the calculated Pearson correlation coefficient, rxy, between T0 and 〈 〉 are also shown in the figures.  (5)) and presented in Table 2. (b)-Dependence of T0 versus 〈 〉 parameters (•) in region 〈 〉 > 0.46 obtained from fits with thermodynamically consistent Tsallis functions with transverse flows (Equation (5)) and presented in Table 2. (c)-The same as in panel (a), but obtained from fits with non-consistent Tsallis functions with transverse flows (Equation (4)) and presented in Table 3. (d)-The same as in panel (b), but obtained from fits with non-consistent Tsallis functions with transverse flows (Equation (4)) and presented in Table 3. The 1-sigma confidence ellipse (corresponding to a 68% confidence interval) of the covariance of parameters T0 and 〈 〉 and the calculated Pearson correlation coefficient, rxy, between T0 and 〈 〉 are also shown in the figures. Now, we can directly compare the results for the threshold border values of 〈 〉, 〈 / 〉, and 〈 〉 estimated in the present analysis in Pb + Pb collisions at = 5.02 TeV, with those extracted in a recent study [24] for Xe + Xe collisions at = 5.44 TeV, where the same theoretical model functions, methods, particle species, and identical intervals were used. The results for the threshold border values of 〈 〉, 〈 / 〉, and 〈 〉 estimated in Pb + Pb and Xe + Xe collisions [24] along with the ratios of the corresponding quantities for two collision types are presented in Table 6. As can be observed in    (5)) and presented in Table 2. (b)-Dependence of T 0 versus β T parameters (•) in region β T > 0.46 obtained from fits with thermodynamically consistent Tsallis functions with transverse flows (Equation (5)) and presented in Table 2. (c)-The same as in panel (a), but obtained from fits with non-consistent Tsallis functions with transverse flows (Equation (4)) and presented in Table 3. (d)-The same as in panel (b), but obtained from fits with non-consistent Tsallis functions with transverse flows (Equation (4)) and presented in Table 3. The 1-sigma confidence ellipse (corresponding to a 68% confidence interval) of the covariance of parameters T 0 and β T and the calculated Pearson correlation coefficient, r xy , between T 0 and β T are also shown in the figures. Now, we can directly compare the results for the threshold border values of N part , dN ch /dη , and β T estimated in the present analysis in Pb + Pb collisions at √ s nn = 5.02 TeV, with those extracted in a recent study [24] for Xe + Xe collisions at √ s nn = 5.44 TeV, where the same theoretical model functions, methods, particle species, and identical p T intervals were used. The results for the threshold border values of N part , dN ch /dη , and β T estimated in Pb + Pb and Xe + Xe collisions [24] along with the ratios of the corresponding quantities for two collision types are presented in Table 6. As can be observed in Table 6  This result is quite interesting, and can be understood if one recalls that the volume of a nucleus is directly proportional to the number of nucleons, that is, to the mass number, A, of the nucleus. Similarly, the number of participant nucleons at some collision centrality is proportional to the volume of the overlap region of two identical colliding nuclei (such as in Pb + Pb and Xe + Xe collisions), which, in turn, should also be proportional to the volume of the nucleus, that is, to its mass number A.
In the end, it is interesting to mention the method of generic (non)extensive statistics (in which the underlying system autonomously manifests its extensive and non-extensive nature), which has been used to analyze the high-energy collision data at the LHC in Refs. [96][97][98] in comparison to the results of extensive Boltzmann and non-extensive Tsallis statistics. In Ref. [96], the yields and their ratios for the identified particles, measured by ALICE Collaboration, were described well using additive thermal approaches taking into account the missing resonance states, van der Waals repulsive interactions, finite volume corrections, and pion chemical potentials, as well as the light and strange quark occupation factors. In Ref. [97], the generic (non)extensive statistics approach was used to describe the transverse momentum spectra of strange hadrons in Pb + Pb collisions at √ s nn = 2.76 TeV, p + Pb collisions at √ s nn = 5.02 TeV, and in p + p collisions at √ s nn = 7 TeV.
The obtained results were compared to those extracted using non-extensive Tsallis and extensive Boltzmann statistics. The differences in the resulting fit parameters (temperature (T), volume (V), and non-extensive parameter d) among three types of statistics were observed in addition to the significant dependencies of the parameters on the size and type of the collision system [97]. It was determined that non-extensive parameter d decreases with increasing collision centrality and energy, indicating that the system departs from non-extensivity (non-equilibrium) and approaches extensivity (equilibrium) with the Boltzmann statistics becoming the proper statistical approach in describing the system in central collisions. It is important to mention that, in the present study, the non-extensivity parameter q values for all studied particle species, resulting from both consistent and non-consistent Tsallis functions with transverse flows, demonstrated a systematic decrease approaching q ≈ 1 value with increasing collision centrality ( N part ). A similar result was obtained for q versus N part dependencies for the analyzed particle species in Xe+Xe collisions at √ s nn = 5.44 TeV in our recent study [24]. Considering that at q ≈ 1 the nonextensive Tsallis distribution reduces to the exponential (equilibrated) Boltzmann-Gibbs distribution, our results for q versus N part dependence agree well with those for the dependence of non-extensive parameter d on collision centrality in Ref. [97]. In Ref. [98], the extensive Boltzmann-Gibbs (BG) and non-extensive Tsallis statistics, generic (non)extensive statistics, and other empirical models were used for the description of the transverse momentum spectra of the charged pions, kaons, and protons produced in p + p and A + A collisions in a wide interval of collision energies. The analytical expressions for the dependencies of the fit parameters on the energy, size of a collision system, and type of the used statistics were obtained [98]. It was deduced that the Tsallis distribution function is more effective within the lower p T range compared to the higher p T interval, and that the Tsallis statistics are more suited for p + p than A + A collisions, while the Boltzmann statistics better describes A+A collisions compared to p + p interactions [98]. It was concluded [98] that the generic (non)extensive statistics was well suited for describing the p T spectra of the charged particles and their antiparticles in both p+p and A+A collisions in a wide range of collision energies. As compared to the pure Tsallis statistics analysis (without embedded transverse flow) used for the comparison with other models in Ref. [98], for the description of p T spectra of the identified particles at different centralities of Pb + Pb collisions at √ s nn = 5.02 TeV in the present study, we used both thermodynamically non-consistent as well as thermodynamically consistent Tsallis distribution functions with incorporated transverse flows (presented in Equations (4) and (5)), which were derived and applied to successfully reproduce the transverse momentum spectra of identified particles in a range up to p T = 5 GeV/c at various centralities of Xe + Xe collisions at √ s nn = 5.44 TeV in our recent study [24]. As observed in Figures 1 and 2, both thermodynamically non-consistent as well as thermodynamically consistent Tsallis functions with embedded transverse flows could reproduce quite well the p T spectra of the charged pions, kaons, and (anti)protons in the region up to p T = 5 GeV/c at various centralities of Pb + Pb collisions at √ s nn = 5.02 TeV in the present study.

Summary and Conclusions
We analyzed the mid-y p T distributions of the charged pions, kaons, and protons and antiprotons measured by ALICE Collaboration at ten centrality classes of Pb + Pb collisions at √ s nn = 5.02 TeV, applying combined minimum χ 2 fits with the non-consistent as well as thermodynamically consistent Tsallis functions with transverse flow. Combined fits of mid-y p T distributions of (π + + π − ), (K + + K − ) and (p + p) with both non-consistent and thermodynamically consistent Tsallis functions with transverse flows reproduced quite satisfactorily the spectra of analyzed particle species in all studied groups of Pb + Pb collision centrality. Combined fits by both consistent and non-consistent Tsallis functions with transverse flow resulted in similar collision centrality ( N part ) dependencies of the parameters T 0 , β T , and q, obtained in Pb + Pb collisions at √ s nn = 5.02 TeV. The parameter β T demonstrated a fairly large growth rate at N part < 71 ± 7 and substantially lower rate of increase at N part > 71 ± 7. This observation was verified quantitatively: the single-power function failed to describe β T versus N part dependence in the whole N part range, whereas the two-power function with two different exponent parameters in regions N part < 71 ± 7 and N part > 71 ± 7 reproduced this dependence very well.
The significantly differing growth rates of β T in regions N part < 71 ± 7 ( dN ch /dη < 251 ± 20) and N part > 71 ± 7 ( dN ch /dη > 251 ± 20) with T 0 staying constant within uncertainties in region N part > 71 ± 7 ( dN ch /dη > 251 ± 20) probably indicate that N part ≈ 71 ± 7 ( dN ch /dη ≈ 251 ± 20) is a threshold border value for a crossover transition from a dense hadronic state to the QGP phase (or mixed phase of QGP and hadrons) in Pb + Pb collisions at √ s nn = 5.02 TeV. The threshold border value β T ≈ 0.46 ± 0.03 (corresponding to N part ≈ 71 ± 7 and dN ch /dη ≈ 251 ± 20), estimated in Pb + Pb collisions at √ s nn = 5.02 TeV in the present analysis, agreed well with the corresponding border value β T ≈ 0.44 ± 0.02, recently obtained in Xe + Xe collisions at √ s nn = 5.44 TeV in Ref. [24], and with almost constant β T values extracted [20] in the BES program at the RHIC in central Au + Au collisions in √ s nn = 7.7 − 39 GeV energy range, where the threshold for QGP production was reached [20,[90][91][92][93][94]. A(132 Xe ) ≈ 1.6 was satisfied. The characters of T 0 versus β T dependence were found to be considerably different in regions β T < 0.46 and β T > 0. 46. The negative correlation between T 0 and β T was quite strong in region β T < 0.46. In contrast to the strong negative correlation in region β T < 0.46, the correlation between parameters T 0 and β T was significantly positive in region β T > 0.46. Hence, it was observed that not only does parameter β T show significantly differing growth rates in regions N part < 71 ± 7 ( dN ch dη < 251 ± 20) and N part > 71 ± 7 ( dN ch dη > 251 ± 20), but also the character of correlation between parameters T 0 and β T differs considerably in the corresponding intervals β T < 0.46 and β T > 0.46. This finding further supports our conjecture that N part ≈ 71 ± 7 ( dN ch /dη ≈ 251 ± 20) is probably a threshold border value for a crossover transition from a dense hadronic state to the QGP phase (or mixed phase of QGP and hadrons) in Pb + Pb collisions at √ s nn = 5.02 TeV.
The non-extensivity parameter q demonstrates a systematic decrease with increasing N part (collision centrality) for all studied particle species. This observation is in agreement with the known fact that the level of thermalization of the system increases with the increasing centrality of heavy-ion collisions at high energies. On the whole, the gap between q(mesons) and q(baryons) decreases with an increase in Pb + Pb collision centrality, and q(mesons) practically coincides within uncertainties with q(baryons) in central collisions with N part > 160. This could indicate quite a large degree of equilibrium and thermalization of QGP produced in central Pb + Pb collisions at √ s nn = 5.02 TeV with N part > 160.