Probing chemical freeze-out criteria in relativistic nuclear collisions with coarse grained transport simulations

We introduce a novel approach based on elastic and inelastic scattering rates to extract the hyper-surface of the chemical freeze-out from a hadronic transport model in the energy range from E$_\mathrm{lab}=1.23$ AGeV to $\sqrt{s_\mathrm{NN}}=62.4$ GeV. For this study, the Ultra-relativistic Quantum Molecular Dynamics (UrQMD) model combined with a coarse-graining method is employed. The chemical freeze-out distribution is reconstructed from the pions through several decay and re-formation chains involving resonances and taking into account inelastic, pseudo-elastic and string excitation reactions. The extracted average temperature and baryon chemical potential are then compared to statistical model analysis. Finally we investigate various freeze-out criteria suggested in the literature. We confirm within this microscopic dynamical simulation, that the chemical freeze-out at all energies coincides with $\langle E\rangle/\langle N\rangle\approx1$ GeV, while other criteria, like $s/T^3=7$ and $n_\mathrm{B}+n_\mathrm{\bar{B}}\approx0.12$ fm$^{-3}$ are limited to higher collision energies.


Introduction
The collision of heavy ions in today's largest particle accelerators provides an excellent tool to explore nuclear and sub-nuclear matter under extreme conditions as they occur e.g. in neutron stars, around black holes or in the early universe. Matter created under these conditions sustains tremendous temperatures, pressures and densities in volumes on the order of ∼ 1000 fm 3 over timescales of 10 −23 s.
Since separating quarks creates new quark-anti-quark pairs from the vacuum in order to bind to color neutral hadrons, a direct measurement of the inner degrees of freedom of strongly interacting matter in heavy ion collisions is, unfortunately, still impossible. In contrast to QED, this phenomenon called confinement forbids perturbative calculations at small momenta. Access to the early and intermediate stage of a collision can be gained via electromagnetic probes e.g. with real photons and virtual photons in the di-lepton channel [1,2,3], via the study of hadrons produced at the chemical freeze-out [4,5], or by flow observables, like v 1 , v 2 , v 3 , · · · .
Typically, the final state of a heavy ion collision involves two parts: 1.) the chemical freeze-out and 2.) the kinetic (or thermal) freeze-out. At the chemical freeze-out inelastic flavour changing reactions cease (particle yields get fixed) and at the kinetic freeze-out also elastic reactions cease (particle 4-momenta get fixed) and the system decouples. This behavior is also reflected in the scattering rates as shown in [6] where the chemical freeze-out is estimated to occur at τ chem = 6 fm. It has been shown, that the kinetic freeze-out proceeds in a more continuous fashion [7] than being an instantaneous process. Assuming thermal equilibrium, particle spectra or ratios can be fitted with statistical and blastwave models to extract kinetic or chemical freeze-out properties like the temperature T and the baryo-chemical potential µ B . Thermal fits of particle ratios resulting from CERN/SPS [8,9,10,11], BNL/AGS [12,13,14] and GSI/SIS [15,16] measurements led to a unified description of the chemical freeze-out line via a hadronic gas model at a constant energy per particle of E / N = 1 GeV [17]. Besides this, other quantities were proposed to define the chemical freeze-out, e.g. a constant total baryon density n B + nB = 0.12 fm −3 [18] or a constant entropy per temperature s/T 3 = 7 [19].
The present investigation introduces a novel approach to determine the chemical freeze-out hyper-surface directly from a combined UrQMD/coarse-grained framework. The full time evolution of Au+Au collisions in the GSI/FAIR and RHIC BES-II energy regime is analyzed for this pur-pose and every final state pion is traced back to the spacetime point of its creation taking into account absorption and decay processes as e.g. N + π ↔ ∆ which affect the freeze-out coordinates. These space-time points define the hyper-surface of the chemical freeze-out of the pions. We focus on pions, because they are the most abundant hadrons and they are produced in sufficient amount at all investigated energies. To obtain the thermodynamic parameters (T , µ B ) on this hyper-surface, the UrQMD data is coarse grained and supplemented by a Hadron Resonance Gas EoS [20]. The energy dependence of the calculated thermodynamic variables is then investigated to test different suggested criteria for the chemical freeze-out.

The model
The present study uses the Ultra-relativistic Quantum Molecular Dynamics (UrQMD) [21,22] transport model in cascade mode. UrQMD is used to compute both the bulk evolution of the system and to determine the space-time coordinates of the chemical freeze-out of the hadrons. We recall that UrQMD is a hadron cascade model that simulates the dynamical evolution of heavy ion collision events by following the propagation of the individual hadrons, modeling their interactions via the excitation of color fluxtubes (strings) and by further elastic and inelastic scatterings. A transition to a deconfined stage is not explicitly included in the cascade mode employed here.

The coarse graining approach
The UrQMD coarse-graining approach [7,23,24,25,26,27] consists in computing the temperature and the baryon chemical potential from the average energy-momentum tensor and net baryon current of the hadrons formed in a large set of heavy ion collision events with the same collision energy and centrality. The computation is done in the cells of a fixed spatial grid at constant intervals of time. In the present study, the cells are four-cubes with spatial sides of length ∆x = ∆y = ∆z = 1 fm and ∆t = 0.25 fm length in time direction. First, we evaluate the net-baryon four current j µ B as and the energy momentum tensor T µν as in which ∆V is the volume of the cell, B i is the baryon number and p µ i stands for the µ component of the four momentum of the hadron i. The sums run over all hadrons N h in the cell. We adopt the Eckart's frame definition [28] and we obtain the fluid four velocity u µ from j µ B as in which γ is the Lorentz factor and v the fluid velocity in natural units (c = = 1). By a Lorentz transformation of the net-baryon current and of the energy momentum tensor in the Local Rest Frame (LRF) of the fluid, we compute the baryon density ρ B and the energy density ε as: Often the chemical freeze-out occurs in cells with an anisotropy between the pressure in the parallel (P ) and in the transverse direction (P ⊥ ), with respect to the beam axis. To take into account this condition, we rescale ε [29,30,24,31] as: where χ = (P ⊥ /P ) 4/3 and This corrections was also applied in all our previous studies [24,25,26,27]. The final step in the coarse graining procedure consists in associating to each cell of the coarse grained grid the temperature T (ε corr , ρ B ) and the baryon chemical potential µ B (ε corr , ρ B ) through the interpolation of a tabulated Hadron Resonance Gas EoS [20].

Chemical freeze-out in the UrQMD model
To test this novel way of tracking down the chemical freezeout coordinates in a full transport simulation we focus on pions. The reason for this is twofold: a) pions are very abundant hadrons in the investigated energy regime and b) pions are stable particles under strong interactions, both facts simplify the reconstruction of the chemical freezeout coordinates with high accuracy. When does a π freezeout chemically? Typically, pions are either produced directly in a string decay (dominant at higher energies) or via N + N → N + ∆, and subsequently ∆ → N + π reactions. Of course the ∆ can be replaced by other resonances. In addition cascades like ∆ → N + ρ, and subsequently ρ → ππ are possible. Not all pions produced initially make it to the final state due to absorption processes, e.g. π + N → N * → K + Λ. To extract the spacetime point of the production of a finally observed π, we follow all observed pions backwards through the evolution until we reach their point of production (i.e. the point at which either their mother-resonances were formed or the pions were produced directly e.g. from a string). This defines the chemical freeze-out coordinates (t, r) for each individual pion.

Results
The results are obtained by analyzing central Au+Au collisions calculated with the UrQMD model. The simulations are used as input to extract the chemical freeze-out coordinates and calculate the fields (T (t, r), µ B (t, r)) with the coarse-graining procedure. Central events are selected via an impact parameter cut at b max = 3.4 fm at all investigated energies. We focus on the chemical freeze-out of pions which finally decouple in the volume with |z| ≤ 5 fm.

Chemical freeze-out times
Let us start with the distribution of the chemical freezeout times at a specific energy to first analyze the influence of the different reconstruction scenarios. For this purpose we show the contributions to the full chemical freeze-out distribution at √ s NN = 19.6 GeV in Fig. 1. The total reconstructed π distribution is shown as a solid black line. It consists of pions that are created hidden in carriers (e.g. ∆, ρ, etc.), shown as red lines and those created directly as pions, shown as blue line. In general the full distribution shows two features: one local maximum centered at ≈ 4 fm and a second small bump arising at ≈ 10 fm. The first peak consists mainly of hadrons created in string excitation reactions, while the second bump consists of decays of resonances and secondary strings. To better interpret  The curves corresponding to the RHIC BES-II energies generally show the same structure: First, a strongly emphasized peak centered at ≈ 5-10 fm (depending on the energy) and second a small bump which shifts towards later times with increasing energy, which is however not visible in the GSI/FAIR energy regime. Here, only the pronounced maximum appears. It is interesting that the peak of the body of the distribution does only very weakly vary with the collision energy in the explored energy regime. Pions created here will mainly evolve through several decay and re-formation cycles of short lived resonances like the ∆ or the ρ. When comparing our chemical freeze-out time distribution with the one given in [7] where the authors investigated the kinetic freeze-out, we clearly observe that the chemical freezeout exhibits a more narrow distribution than the kinetic freeze-out. This finding validates therefore that the conceptual idea of an instantaneous chemical freeze-out is to some extent reasonable while the kinetic freeze-out studied in e.g. [7] is shown to be a more continuous and spread-out process. The position of the small bump in contrast depends on the energy. We relate this effect to the γ factor of resonances formed in string decays. If e.g. a ∆ resonance is formed from a di-quark within a string, which decays into a ρ + N and subsequently into a final state pion, the creation time satisfies: t chem. ≈ √ s NN /(2m p ) + Γ −1 ∆ . The average chemical freeze-out times τ chem can be found in Tab. 1 stating that the chemical freeze-out in the RHIC BES-II energy regime occurs at ≈ 7 fm.

Temperature and chemical potential distributions
Now that the freeze-out four-coordinates are determined, the question arises how the coarse-grained temperatures and the baryon chemical potentials will be distributed. For this we show in Fig. 3 the distribution of the coarsegrained chemical freeze-out temperatures of pions with |z| ≤ 5 fm in Au+Au collisions from E lab = 1.23 AGeV to √ s NN = 62.4 GeV (same colors as in Fig. 2) from UrQMD as well as we show the distributions of the baryon chemical potential at the chemical freeze-out in Fig. 4. We observe two major features in the distributions of temperatures extracted from the reconstructed hypersurface: a) with increasing beam energy, the temperature distributions merge into one peak at ≈ 150 MeV (in line with expectations from a thermal model analysis) while the typical widths of the temperature spread is FWHM ≈ 50 MeV, different from the expectations from a thermal model analysis with a single temperature. b) The distributions of the baryon chemical potentials at chemical freeze-out is very narrow at each energy. As expected the chemical potential decreases with increasing energy. In contrast to the temperatures we do not observe a saturation of µ B , but a continuous decrease.

Energy dependence
In this section the collision energy dependence of the temperature and the baryon chemical potential are investigated on the chemical freeze-out hyper surface. We compare the chemical freeze-out values to kinetic freeze-out values. The kinetic freeze-out is defined as the point of the particles last interaction and can also easily be extracted from the presented algorithm. We show in Fig. 5 the average temperature at the kinetic freeze-out (blue circles) and the average temperature at the chemical freeze-out (red circles) as a function of the collision energy. Firstly, it is noticeable and important that T chem > T kin over the whole energy range. This implies that indeed the chemical and the kinetic freeze-out processes happen sequentially. The energy dependence of the absolute difference ∆T = T chem − T kin between the two freeze-out temperatures stays according to the present analysis approximately constant (20 ± 5 MeV). A saturation of the chemical freeze-out temperature is observed at ≈ 150 MeV while the kinetic freeze-out temperature saturates at ≈ 130 MeV. Fig. 6 shows the average baryon chemical potentials at the chemical and the kinetic freeze-out. Both quantities drop rapidly with increasing energy, lining up with the vanishing chemical potential inferred from observed baryon-charge symmetric matter formed at top RHIC and LHC energies.

Phase diagram
After the analysis of the time distributions and the energy dependence of the chemical freeze-out quantities (temper-atures and baryon chemical potentials), we will now relate our results to the phase structure of QCD in the phase diagram of nuclear matter. Fig. 7 shows the calculated average temperature T and the average chemical potential µ B as red circles for the chemical freeze-out and as blue circles for the kinetic freeze-out. Fits to the experimental data points obtained at RHIC [32,33,34,35,36], CERN/SPS [37,38], BNL/AGS [37,38], GSI/SIS [39,40,41] and the recent HADES point [42] are shown as green stars (summarized in [43,44]). We observe in general that our calculated chemical freeze-out line matches the estimates from thermal model fits to the experimental data surprisingly well. This result is remarkable because it is ex-ante not expected that a purely hadronic transport model, which does neither involve the concept of the chemical freeze-out explicitly nor a phase-transition to a deconfined phase, should reproduce (T , µ B ) combinations near the data points obtained from thermal model fits. This raises a question: Why does a non-equilibrium transport model without phase-transition and chemical break-up reproduce thermal (i.e. equilibrium) quantities such precise? A hint towards the answer can be found the scattering rates. It is well known from e.g. chemistry that the onset of equilibrium can be characterized by the ratio of the expansion rate to the scattering rate. This means to maintain equilibrium, the (inelastic) scattering rate Γ must be much larger than the expansion rate θ. Where θ/Γ = Kn, can be interpreted as the Knudsen number Kn. Typically, Γ ∼ f i f j σ ij , with f i being the phase space density of species i and σ ij being the (inelastic) interaction cross section, and θ = ∂ µ u µ , the divergence of the 4-velocity field [45]. I.e., if the scattering rate is larger than the expansion rate then enough equilibrating reactions happen to keep the system in equilibrium. If otherwise the scattering rate is smaller than the expansion rate then the system does not have enough time to adjust its properties and the state becomes frozen for the remaining time of the expansion. While in our analysis, the competition between scattering rate and expansion rate is modelled by microscopic dynamics, it can also be used to calculate the chemical and kinetic freeze-out in a self consistent way in expanding gases and fluids, see e.g. [45,46]. Therefore, the temperatures and the baryon chemical potentials at chemical and kinetic freeze-out as calculated here emerge from the local interplay of the elastic and inelastic collision rates with the local expansion rates of the system. Thus, the chemical and kinetic freeze-out lines are not related to a phase-transition or a cross over.

Testing chemical freeze-out criteria
Our analysis can be further deepened by the investigation of quantities which have been proposed to classify the chemical freeze-out line. The systematic investigation of the chemical freeze-out in experiments and theory led to the proposal of several quantities which apparently seem to function as criteria for the chemical freeze-out to occur: . Also shown are thermal model fits to RHIC [32,33,34,35,36], SPS [37,38], AGS [37,38], SIS data [39,40,41] and the recent HADES point [42] as green stars.
[17] the average energy per particle E / N was proposed to stay constant along the chemical freezeout line with E / N = 1 GeV. 2. Secondly, an entropy based criterion, s/T 3 = 7 [19] was suggested for the meson dominated energy regime. 3. It was further proposed that the total baryon number can also serve as an estimation for the freeze-out line by demanding n B + nB = 0.12 fm −3 [18].
These suggestions can be directly accessed and tested in the coarse-graining method employed here.

Average energy per particle
With increasing collision energy a major part of the additional energy is used to produce new (heavier) particles than to increase the momentum of the existing particles [47] leading to the limiting Hagedorn temperature for hadrons. In fact, the study of hadronic abundances at SIS, AGS, SPS, RHIC and LHC energies has shown that the chemical freeze-out line can be characterized by a constant E / N of ≈ 1 GeV [17] known as the Cleymans-Redlich criterion. This pioneering observation induced many other works to investigate measured hadron ratios which concluded to similar values ranging between E / N = 0.96 -1.08 GeV [43,48]. In [38] using a statistical model with two thermal sources (TSM) it was predicted that the average energy per particle as a function of √ s NN may not be constant, but depend on the resonance share of the produced matter and suggested an increase of E / N to 1.1 -1.2 GeV for √ s NN ≤ 10 GeV.
We show in Fig. 8 the calculated E / N ratio on the chemical and kinetic freeze-out surface arising from the coarse-graining method in dependence of the collision energy. The results at the chemical freeze-out are shown as red circles and the results for the kinetic freeze-out are shown as blue circles. First of all, the computed average energy per particle at the chemical freeze-out hyper surface varies between 0.95 and 1.2 and is therefore very close to the phenomenological 1 GeV/particle. In line with [38,48], we observe a slight energy dependence of the ratio with an increase towards 1.2 at the baryon dominated energy region. When studying this ratio on the kinetic freeze-out hyper surface, it is surprising that the same quantity computed here yields only a slightly lower value of E / N ≈ 0.95 GeV. This can be explained when examining how E / N varies with time which has been done e.g. in [6]. The system is characterized by a plateau in the energy per particle over the duration of the chemical equilibration. Thus we can conclude that the chemical freeze-out line is indeed characterized by E / N ≈ 1 GeV/particle. Towards lower energies, we predict an increase by about 20 % in line with [38,48].

Entropy density
In thermodynamics, the entropy density is the determining quantity that characterizes when (chemical) equilibrium sets in and is thus a prime candidate to constrain the chemical freeze-out. For an ideal gas of massless hadrons with zero net-baryon density the dimensionless term s/T 3 is equivalent to the number of active degrees of freedom. E.g. a pion gas with s/T 3 = 3 would have exactly 3 degrees of freedom. In reality, however, pions are not the only degrees of freedom in the hadron gas, and the entropy has contributions from all mesons and baryons, thus s/T 3 is a proxy for the effective degrees of freedom. This relation was extensively studied with the conclusion that chemical equilibration is characterized by a constant entropy density per cubic temperature of 7 (equivalent to 7 effective degrees of freedom) by thermal model studies [19]. Therein it was further shown that the decomposition of s/T 3 = 7 into hadronic and mesonic contributions reveals an interchange in the dominant contribution around √ s NN = 7.7 GeV. At low energies the degrees of freedom are baryon dominated while at high energies mesonic degrees of freedom are the dominating player. The calculated entropy density s/T 3 is shown in Fig.  9. The red circles show the values obtained in the present analysis at the chemical freeze-out and the blue circles correspond to the values at the kinetic freeze-out, both in dependence of the collision energy. The filled circles represent the meson dominated region and the empty circles show the baryon dominated region. To avoid numerical problems related to the temperature extraction in cells with high baryon density and low energy density (i.e. cold cells), we include only cells with T ≥ 30 MeV in the averaging. Extracting the chemical freeze-out from UrQMD, we find that between 7.7 GeV and 62.4 GeV s/T 3 remains at values between 6 and 7 lining up with suggested values of Ref. [19]. Below the √ s NN = 8 GeV, the entropy per cubic temperature rises rapidly as expected from [19]. For the kinetic freeze-out, we find that s/T 3 saturates be- tween 4 − 5 in the meson dominated region and rises in the baryon dominated region. In conclusion, the s/T 3 criterion motivated by the number of degrees of freedom is found to be in line with the previously obtained value of 7 [19] in the whole energy range and above √ s NN = 8 GeV.

Baryon and anti-baryon density
Finally, we explore the suggested baryon-anti-baryon criterion with coarse grained UrQMD. Fig. 10 shows the calculated n B +nB in dependence of the collision energy. The results for the chemical freeze-out are shown as red circles, while the results for the kinetic freeze-out are shown as blue circles. The results for the total baryon and antibaryon density at the chemical freeze-out start at a density of 0.35 fm −3 at the low energies and drop rapidly to 0.15 fm −3 at √ s NN = 20 GeV. This increase at low energies is qualitatively in line with expectations in [43]. The results at the kinetic freeze-out follow the same trend, but the decrease is more moderate starting at 0.15 fm −3 and ending at 0.05 fm −3 . Generally, the baryon density criterion shows a stronger energy dependence than the previously discussed criteria for the chemical freeze-out line. Nevertheless, towards higher energies ( √ s NN ≥ 20 GeV) also this criterion can be applied to characterize the chemical freeze-out line.

Conclusion
In this article we have developed a novel approach to determine the chemical freeze-out hyper-surface directly from a microscopic simulation. The UrQMD transport model was employed to simulate the underlying events and the microscopic evolution of the system. Using a new algorithm, we traced final state particles back to their original creation space-time coordinate through sequential decay and re-formation chains allowing to reconstruct the chemical freeze-out hyper-surface using a coarse-graining method. We found that the average chemical break-up time remains constant at ≈ 7 fm above √ s NN = 7.7 GeV.
The extracted temperature and baryon chemical potential values follow the trend of the established thermal model fits to experimental data. Since UrQMD neither involves a phase-transition to a deconfined state of matter nor the explicit concept of the chemical freeze-out, our results indicate that the chemical freeze-out line may not be indicative of the QCD phase transition, but is defined by the competition of the inelastic scattering rate with the expansion rate. The investigation of previously suggested freeze-out criteria reveals that indeed a constant s/T 3 = 7 and n B + nB ≈ 0.12 fm −3 are good proxies to characterize the chemical freeze-out curve at high collision energies ( √ s NN ≥ 20 GeV). The originally suggested Cleymans-Redlich criterion E / N = 1 GeV, however, provides the best estimate for the chemical freeze-out line over all energies.