Hadronization conditions in relativistic nuclear collisions and the QCD pseudo-critical line

We compare the reconstructed hadronization conditions in relativistic nuclear collisions in the nucleon-nucleon centre-of-mass energy range 4.7-2760 GeV in terms of temperature and baryon-chemical potential with lattice QCD calculations, by using hadronic multiplicities. We obtain hadronization temperatures and baryon chemical potentials with a fit to measured multiplicities by correcting for the effect of post-hadronization rescattering. The post-hadronization modification factors are calculated by means of a coupled hydrodynamical-transport model simulation under the same conditions of approximate isothermal and isochemical decoupling as assumed in the statistical hadronization model fits to the data. The fit quality is considerably better than without rescattering corrections, as already found in previous work. The curvature of the obtained"true"hadronization pseudo-critical line kappa is found to be 0.0048 +- 0.0026, in agreement with lattice QCD estimates; the pseudo-critical temperature at vanishing mu_B is found to be 164.3+-1.8 MeV.


I. INTRODUCTION
It is the goal of Quantum Chromo-Dynamics (QCD) thermodynamics to study the phase diagram of strongly interacting matter. Its most prominent feature, the transition line between hadrons and partons, in the plane spanned by the baryon-chemical potential µ B and the temperature T , is located in the non-perturbative sector of QCD. Here, the theory can be solved on the lattice and has recently led to calculations of the curvature of the parton-hadron boundary line [1][2][3][4][5][6][7][8]. This line can also be studied experimentally, in relativistic collisions of heavy nuclei, where apparently local thermodynamical equilibrium is achieved at a temperature well above the (pseudo-)critical QCD temperature T c . Expansion and cooling then take the system down to the phase boundary where hadronization occurs. We have lately demonstrated [9][10][11] that post-hadronization inelastic rescattering, chiefly baryon-antibaryons annihilation, is an important feature of the process, which drives the system slightly out of equilibrium from the primordial hadronization equilibrium, implying an actual distinction between hadronization and chemical freeze-out. This rescattering stage is taken into account in state-of-the art simulations of the QGP expansion [12][13][14], where the local equilibrium particle distribution (through the so-called Cooper-Frye formula) at some critical values of T and µ B is used to generate hadrons and resonances which subsequently undergo collisions and decay. By calculating the modification of the multiplicities brought about by the rescattering stage -the so-called afterburning -it is possible to reconstruct the hadronization point by means of a fit to the multiplicities in the framework of the Statistical Hadronization Model (SHM). Strictly speaking, this method allows to pin down the latest chemical equilibrium point [11] henceforth denoted as LCEP -i.e. the point when the primordial chemical equilibrium starts being distorted by the afterburning. As equilibrium is an intrinsic feature of hadronization [15][16][17] -as shown by the analysis of elementary collisions -most likely LCEP coincides with hadronization itself, as the maintaining of full chemical equilibrium in a rapidly expanding hadronic system, for the time needed to produce a measurable temperature shift, is highly unlikely. As the primordial system temperature (baryon-chemical potential) shifts upward (downward) with increasing collision energy, an ascending sequence of experimental energies can, thus, map a sequence of LCEPs or hadronization points along the QCD transition line. This was the main point of ref. [10] where we showed that, indeed, the reconstructed LCEPs seem to follow the extrapolated lattice QCD pseudo-critical line in the (µ B , T ) plane as determined in, e.g., ref. [3]. The agreement between lattice QCD calculations and the reconstructed hadronization points in relativistic heavy ion collisions seem to imply that, in the examined energy range ( √ s NN > 7 GeV), the pseudo-critical line has indeed been crossed, and, thus, those energies lie above the so-called "onset of deconfinement" [18]. This conclusion is less straightforward than it may seem at a glance because, as has been mentioned, hadron formation is evidently a universal statistical process [15][16][17] in all kinds of collisions at the same hadronization temperature, with a difference in the strangeness sector, whose phase space appears to be only partially filled in elementary collisions [19][20][21][22][23] 1 . Indeed, even if the strangeness phase space was fully saturated in elementary collisions, if hadron production process in a nuclear collision was fully consistent with a picture of subsequent and independent elementary hadronic reactions, strange particle production would be strongly suppressed by the exact strangeness conservation over the typical small volumes of an elementary collision (canonical suppression) and subsequent hadronic inelastic collisions would not be able to raise multi-strange particle abundance to the measured one, as it is shown by transport calculations [24][25][26]. Hence, the observation of a strange particle production in agreement with the prediction of the SHM for a coherent, large volume, and the agreement between lattice QCD extrapolations of the pseudo-critical line and the reconstructed hadronization point is a strong evidence that the pseudo-critical line has been overcome. However, this must cease to happen at some sufficiently low centre-of-mass energy, implying the failure of at least one of the above conditions. Estimates remain uncertain at present, pointing to a region between 4 and 8 GeV.
In this paper, we reexamine the hadronization conditions in relativistic heavy ion collisions over the energy range from low SPS (7.6 GeV) to LHC (2.76 TeV) by using hadronic multiplicities. We also include the highest AGS energy point at √ s NN = 4.5 GeV to probe the aforementioned deconfinement conditions further down in energy. For this purpose, we take advantage of an improved initialization of the afterburning process by enforcing a particle generation stage -or hydrodynamical decoupling -in UrQMD [27][28][29][30] at fixed values of energy density corresponding to mean temperatures and chemical potentials equals to those determined in ref. [10] and show how this leads to a further remarkable improvement of the fit quality compared to the plain statistical model fits [10,11,31]. Finally, we compare the resulting curvature of the LCEP-hadronization curve in the (µ B , T ) plane with the predictions of lattice QCD, reporting a good agreement.

II. AFTERBURNING AND MODIFICATION FACTORS
As has been mentioned, we studied the effect of post-hadronization rescattering on hadron multiplicities and on the associated SHM fits in previous publications [9][10][11] employing a hybrid model [34,35] with a hydrodynamic expansion of the QCD plasma terminated at a predefined point where local equilibrium particle generation is assumed (hadronization), followed by a hadronic rescattering stage modelled by UrQMD [29] (afterburning). For the fluid dynamical simulation we employed an equation of state which follows from a so-called "combined hadron-quark model". It is based on a chiral hadronic model which provides a satisfactory description of nuclear matter properties. The quark phase is introduced as a PNJL type model. The transition from the hadronic to the quark phase occurs at about T c ≈ 165 MeV for µ B = 0, as a smooth crossover as shown in [32].
For each hadronic species a so-called modification factor is extracted which is defined as the ratio between the final multiplicity after the actual chemical (and kinetic) freeze-out (which is now species-dependent) and its value without afterburning (at hadronization): The modification factors are then used as additional multiplicative factors to the theoretical equilibrium multiplicity yields in the SHM fit, ready to be compared to the data. Note that in the calculation of the modification factors, all weak decays are turned off, but all strong and EM decays are turned on; this limits the data analysis to measurements of multiplicities corrected for the weak decay feed-down. In our previous studies, the hydrodynamic decoupling procedure was inspired by the so-called "inside-outside cascade" mechanism: the transition from the fluid dynamical phase to the hadronic transport part occurs in successive transverse slices, of thickness 0.2 fm, whenever all fluid cells of that slice fall below a critical energy density, that is six times the nuclear ground state density ǫ ≈ 850 MeV/fm 3 . In fact, in the present investigation, we have implemented an approximate isothermal termination of the hydrodynamical stage at some pre-established temperature T CF (the subscript CF stands for Cooper-Frye). This is certainly in much better accordance with the underlying picture of a statistical hadronization as well as with the previously discussed concept of LCEP, which is determined at a fixed value of the proper temperature. For the decoupling, the hypersurface is defined by a fixed energy density -at LHC energy -of approximately 0.360 GeV/fm 3 which corresponds to a mean hadronization temperature close to 165 MeV at zero baryon-chemical potential (for the lower energies, see table I). The cell-to-cell temperature and chemical potential fluctuations corresponding to such a hydrodynamic decoupling procedure, at some given collision energy, are small. For instance, the dispersion of the temperature at LHC energy is ∼ 1.5 MeV, the dispersion of the chemical potential at the SPS energy is of the order of 10 MeV; these values are comparable or smaller than the parameter fit errors (see table V).
The UrQMD model employs the hypersurface finder outlined in ref. [30], which is used in the Cooper-Frye prescription and sampled to produce hadrons in accordance with global conservation of charge strangeness baryon number and the total energy. It should be pointed out that the calculated modification factors do depend on the chosen temperature T CF ending the hydrodynamical expansion [33]. Ideally, this should coincide with the actual T LCEP at each energy, which is a priori unknown, except for a reasonable lower bound set by the chemical freeze-out temperature as determined in the traditional, plain, SHM fit. One may then wonder to what extent T LCEP , which is the outcome of the subsequent SHM fit -corrected for afterburning -is affected by the chosen T CF . In general, larger T CF involve larger deviations of the modification factors from unity, so one could expect that the fit is influenced to such an extent that the corrected SHM fit tends to reproduce the initially chosen value T CF , making the whole method non-predictive. However, this would be the case only if the final particle multiplicities after freeze-out were independent of T CF . It was checked that in hybrid simulations this does not happen and final particle yields do depend on the chosen decoupling condition. In fact it appears that the change in the extracted T LCEP tends to be smaller than in the input T CF (i.e. independent of the exact values of the modification), and T LCEP shows a trend toward a definite value. For instance, for Pb-Pb collisions at √ s NN = 8.7 GeV (see next section for details), for an input T CF = 144 MeV we obtained T LCEP = 155 MeV whereas for T CF = 161 MeV we obtained T LCEP = 163 MeV. In summary, the method converges.
The optimal situation, as has been mentioned, is T CF = T LCEP , which could be achieved by an iterative procedure; however, it would be computationally expensive and not worth the effort when -in view of the above observation -the difference between T CF and T LCEP is only few MeV's. Altogether, the small differences between T CF and T LCEP in our analysis (see table V) in the next section) make us confident that the fitted thermal parameters are fully significant, with only a marginal dependence on the difference (T CF − T LCEP ). The modification factors for the strongly stable hadrons have been calculated according to the method described in the previous section and are quoted in table II . The decoupling conditions in terms of temperature, that is T CF , and baryon-chemical potential are those determined as LCEP's in our previous papers for the most central collisions at the LHC [11] and for SPS points [10] (see table V further below). For the AGS point we did not have any clue about the T LCEP value, so we iterated the procedure until the fitted T LCEP was reasonably close to the T CF . The modification factor for the φ is more difficult to extract than for other, strongly stable, particles as it is not the φ itself which is absorbed, but its decay products which rescatter, making φ reconstruction hardly feasible. Assuming that any rescattering of a decay product will lead to a loss of a φ, the lower bound of the survival probability has been estimated at 2.76 TeV to be about 0.75 [41][42][43]. At all lower energies, for afterburning is expected to be less important for this meson, we have used an educated guess of 0.875 which is the mean value between 0.75 and 1, varying between these bounds in order to check the stability of the best fit solutions.
The particle set used in the analysis is the intersection between the available measured multiplicities and the set of particles for which a modification factor was calculated, see table III. All data refer to central collisions of Au+Au (at the AGS) and Pb+Pb (at SPS and LHC). As has been mentioned, we have confined ourselves to data sets where weak feed-down was subtracted in order to make a proper comparison between corrected (with modification factors) and non-corrected fits.
The SHM, the formulas for primary and final multiplicities, the fitting procedure with and without modification factors have been described in detail elsewhere [9]. Herein, we simply summarize the obtained results in table V. The corrected fits are of remarkable better quality with respect to the plain SHM fits, as it is shown in fig. 1, confirming previous findings. One exception stands out, the SPS point at 8.7 GeV, which is the point where the ratio K + /π + attains its maximum observed value [51]. Indeed, the measured ratio K + /π + overshoots the statistical model a This is the ratio p/π + b Our extrapolation [9] based on measurements in ref. [49] c Our extrapolation [9] of spectra measured in ref. [50]. The NA49 data compilation [44] quotes 4.25±0.28 by M. Utvic. d Ω +Ω e Interpolation to 0-5% centrality quoted in ref. [11] f Number of participants = net baryon number of the fireball prediction [52] by more than 2σ (see table IV), a discrepancy which is not cured by the afterburning correction.
Notably, the χ 2 /dof is of full statistical significance once the modification factors are introduced in the two highest energy points. Moreover, there is a further improvement of the fit quality with respect to the previous, non-isothermal, fits [10,11].
The quoted errors in table V are the fit errors. There are additional small systematic uncertainties on the fit parameters related to the errors on the modification factors. These errors stem from the uncertainties on the crosssections used in UrQMD and from finite Monte Carlo statistics. The former are difficult to estimate, whereas the latter are simpler; in our runs they are of the order of few percent for all particles, thus they do not imply any significant variation of the best fit parameter values. The only largely unknown modification factor is, as has been mentioned, the φ's, for which we could obtain an estimate of 0.75 at the top energy point √ s NN = 2760 GeV. As the rescattering of a neutral meson is expected to diminish in a lower multiplicity environment, one can reasonably set a lower bound of 0.75 at all lower energies. To estimate the effect of the uncertainty, we have varied the φ modification factor to 0.75 and 1 in turn at each energy point. The resulting variation of the fit parameters is within 1 MeV for the temperature and few MeV's for the baryon chemical potentials so that the relative systematic error is always less than 1%, thus below the fit error.

IV. CURVATURE OF THE PSEUDO-CRITICAL LINE
Finally, we have used the LCEP points in table V to fit the curvature of the pseudo-critical line in the (µ B , T ) plane. The curvature parameter κ is defined by the equation: which is the same formula used in lattice calculations. As the LCEP points in (µ B , T ) plane have errors on both coordinates, we have minimized the χ 2 : where In the above equation, µ Bi and T i are the output of the SHM fit for the ith centre-of-mass energy point, while C i is their corresponding covariance matrix. The µ 0 Bi 's are free parameters which represent the "true" values of the chemical potential in the fitted curve, T c (µ 0 Bi ) being the corresponding "true" temperatures. Therefore, the free parameters in this fit are the chemical potentials µ 0 Bi , whose position is strongly constrained by the "measured" values µ Bi , the value of the pseudo-critical temperature T c (0) and κ.
It should be pointed out that the equation (2) is a quadratic approximation of the actual pseudo-critical line, hence deviations are expected at large values of the chemical potentials. Therefore, we have first excluded the lowest energy point and made a fit to the four highest energy points at our disposal. We have also compared with the freezeout points, for which many systematic studies have been done in the past [53]. The fitted values of T c (0) and κ are reported in table VI while the fitted curves are shown in fig. 2. It can be seen that the fit quality is excellent for the LCEP points and satisfactory for the plain freeze-out points. The systematic error on the curvature due to the uncertainty in the φ meson modification factors has been obtained repeating the fit with the varied (µ B , T ) points and turned out to be 0.0006. The lowest energy point falls below the fit curve in both cases. There are three possible explanations for this: 1. the (mundane) effect of having excluded it from the fit; 2. the quadratic approximation in (2) falls short at such large chemical potential values; 3. the lowest energy point did not reach the pseudo-critical transition line, and so the onset of deconfinement can be located between 4.5 and 7.6 GeV.  [5] 0.0149±0.0021 [7] 0.0135±0.0020 [8] 0.020±0.004 The latter hypothesis is indeed the most intriguing, but its very consideration requires the ruling out the first two. If the AGS point is included, the fit quality is not as good as the 4 point fit, still it is within statistical significance, as it can be seen in table VI. Yet, there is some tension between the fitted curve and the two extremal points (LHC and AGS) which both undershoot the curve by 4 and 15 MeV respectively, as it can be seen in fig. 2). On the other hand, including a quartic term λ(µ B /T c (0)) 4 improves the fit but the limited number of points and the limited range does not allow to pin down both the quadratic and the quartic term at the same time; indeed, the fit has multiple solutions and the best fit is awkwardly found for κ ≃ 0 (see table VI). Finally, returning to fig. 2, which is our main result showing the estimated QCD pseudo-critical curve in the (µ B , T ) plane, it is appropriate to compare it with recent lattice QCD calculations (see table VII). Because of the pseudocritical nature of the transition, both T c (0) and κ parameters depend on the observable used to define it [4]. It could be therefore expected that these parameters will somewhat differ from those extracted with the fluctuation of conserved charges [54], which can be directly calculated in lattice QCD but are definitely less robust observables in relativistic heavy ion collisions with respect to mean multiplicities [55].
In our comparison, we have quoted all recent literature on the subject. We find that our main value of 0.0048 is in slightly better agreement with lower estimates [1,2,4]. We also note that our value is compatible with a recent estimate based on a comparison between lattice QCD and data [56].

V. CONCLUSIONS
In summary, we have determined the hadronization conditions (strictly speaking the latest chemical equilibrium points) in relativistic heavy ion collisions, by using an improved calculation of the post-hadronization rescattering correction. The quality of the statistical model fit considerably improves with respect to traditional fits without afterburning corrections as well as with respect to our previous calculations. The pseudo-critical temperature at zero µ B , determined with hadronic multiplicities, turns out to be T c = 164 MeV, which is significantly higher than lattice QCD calculations based on different observables. We find good agreement between the extracted curvature of the hadronization curve and the corresponding QCD lattice calculations for the pseudo-critical line, with a preference for the lower estimates.
At this stage, it is not possible to make a definite statement about the crossing of the pseudo-critical line at the lowest energy point at √ s NN = 4.7 GeV. This is due to lack of appropriate data and this issue will be tackled by future experiment in that energy range (NA61 at SPS and the facilities NICA and FAIR). Our observations might have interesting implications for the location of the critical point [57]; we note that a recent analysis [58] locate it at T ≃ 165 MeV and µ B ≃ 95 MeV which just sits on our hadronization curve.