Ambiguities in the hadro-chemical freeze-out of Au+Au collisions at SIS18 energies and how to resolve them

The thermal fit to preliminary HADES data of Au+Au collisions at $\sqrt{s_{_{NN}}}=2.4$ GeV shows two degenerate solutions at $T\approx50$ MeV and $T\approx70$ MeV. The analysis of the same particle yields in a transport simulation of the UrQMD model yields the same features, i.e. two distinct temperatures for the chemical freeze-out. While both solutions yield the same number of hadrons after resonance decays, the feeddown contribution is very different for both cases. This highlights that two systems with different chemical composition can yield the same multiplicities after resonance decays. The nature of these two minima is further investigated by studying the time-dependent particle yields and extracted thermodynamic properties of the UrQMD model. It is confirmed, that the evolution of the high temperature solution resembles cooling and expansion of a hot and dense fireball. The low temperature solution displays an unphysical evolution: heating and compression of matter with a decrease of entropy. These results imply that the thermal model analysis of systems produced in low energy nuclear collisions is ambiguous but can be interpreted by taking also the time evolution and resonance contributions into account.


I. INTRODUCTION
A recent statistical model analysis [1] of the new HADES collaboration data yields an unexpectedly low temperature best χ 2 fit for hadron yields in a standard chemical freeze-out model. Statistical models are a well established tool to describe hadron production in heavy ion collisions in experiments at the Bevalac, SIS, AGS, SPS, RHIC, and LHC accelerators. It is of great interest that static statistical models with only a handful of thermodynamic parameters provide a surprisingly good description of the particle yields from system with complicated non-equilibrium dynamics and interactions. The ideal hadron resonance gas model (HRG) gives a generally good description of the many experimentally observed hadron yields measured at various collision energies [2][3][4][5][6]. This approach assumes that the chemical composition of the system, and thus the final particle multiplicities, are fixed at late stage of heavy ion collisions, the so-called chemical freeze-out. The chemical freezeout parameters for each collision system are obtained through a fit to measured particle multiplicities. Collision energy dependence of the extracted parameters like temperature and chemical potentials defines the chemical freeze-out line, mapping heavy-ion collision experiments to the QCD phase diagram [7,8]. Although the statistical model has different implementations, which can mainly vary in the way hadronic interactions are treated [9][10][11][12][13], with few exceptions [14], they give very similar chemical freeze-out curves. Consequently, the existence of a universal description of the last point of chemical equilib-rium in heavy ion collisions is of great importance to the interpretation of heavy ion experimental data. Any deviations or discrepancies with this picture would therefore be of great interest, having strong implications on the applicability of the standard paradigm governing high energy heavy ion collisions. This paper shows that the hadron yields in the few GeV collision energy regime pose a conundrum for the thermal model fit. The particle yields measured by the HADES collaboration at SIS18 accelerator at GSI can be described similarly well by two distinct sets of thermodynamic parameters. The same degeneracy is observed in transport simulations within the UrQMD model. Studying the time evolution of the system within UrQMD allows to elaborate on the viability of the two solutions.

II. DATASET USED
The HADES collaboration at the SIS18 accelerator at GIS has measured a set of preliminary hadron multiplicities at 10% most central Au+Au collisions at √ s N N = 2.4 GeV [15][16][17][18]. These data are summarized in Table I. A thermal model analysis of these preliminary HADES data I has been previously published in Ref. [1]. The authors extracted thermodynamic properties of the system created in the heavy ion collisions by assuming that all hadron multiplicities are fixed at a single chemical freezeout. The freeze-out of the system was described within an ideal HRG model. Additionally, it was assumed that  Table I. Different locations of the freeze-out correspond to different considered freezeout scenarios, see text for details. (b): χ 2 profiles of the performed fits. The "no clusters" scenario gives the best fit, however produces two distinct minima. To distinguish between the two minima the π feeddown fraction is presented which is calculated as a ratio of the number of charged pions that stem from resonance decays to the total number of charged pions, (π + + π − ) feeddown /(π + + π − ) total .
light nuclei are not present at the chemical freeze-out stage, but that they are formed instead at a later stage via a different mechanism. The protons calculated at the chemical freeze-out do therefore include both, the protons measured as free protons and those which are later bound into light nuclei. The analysis used p, π +,− , K +,− , and Λ yields as input to the fit. The measured proton yield was calculated as a sum of the directly measured unbound protons and protons bound inside the measured light nuclei. Here the same dataset is used as in Ref. [1].

III. THERMAL MODEL ANALYSIS
The thermal model analysis performed in Ref. [1] suggested that the chemical freeze-out in Au+Au collisions at √ s N N = 2.4 GeV occurs at rather cold and dilute system, characterized by the following set of thermodynamical parameters: T = 49.6 ± 1 MeV, µ B = 776 ± 3 MeV, µ I3 = −14.1 ± 0.2 MeV, µ S = 123.4 ± 2 MeV, and γ s = 0.16± 0.02. In the following first a similar analysis using the same data is performed. Here it is assumed that the particle yields are fixed at the chemical freezeout which is described by ideal HRG in the grand canonical ensemble. The freeze-out state is then related to the thermodynamic parameters of a hadron resonance gas: temperature T , baryon, electric, and strangeness chemical potentials µ B , µ Q , and µ S , respectively. The system size is defined through the freeze-out radius R or volume V = 4 3 πR 3 . Strangeness undersaturation is taken into account via a strangeness suppression factor γ S [19]. The electric and strangeness chemical potentials µ Q and µ S are fixed by the charge (electric and strange) content of the colliding nuclei, namely the total strangeness vanishes n S = 0 and the electric to baryon charge ratio equals to n Q /n B = 0.4. It should be noted that in the Ref. [1] µ Q and µ S were used as free fit parameters, this led to N df = 0 and resulted in a slightly different set of thermodynamical parameters, as compared to the analysis presented here. The ideal HRG model is used, i.e. at the freeze-out the system is represented by a multi-component gas of free hadrons and resonances in equilibrium. The particle list consists of hadrons listed in Particle Data Tables 2020 [20] with an established status. Our analysis is performed using the open source Thermal-FIST package version 1.3 [21]. In an alternative scenario one can assume that the yields of light nuclei are fixed at the chemical freezeout together with all other hadrons. The thermal model works remarkably well for describing the yields of light nuclei across a broad range of collision energies [22][23][24]. Thus, in both the discussed scenarios, the formation of nuclei at a late stage due to coalescence from primordial nucleons as well the direct creation of nuclei at the chemical freeze-out can give similar results on the yields of nuclei. It is challenging to explain how the nuclei may survive the temperatures which are an order higher than their binding energies if the second scenario is the correct one. However, some progress on the understanding of this phenomenon has recently been made [25][26][27]. See also [28][29][30] for possible ways to distinguish between thermal model and coalescence.
To take into account the different possibilities for the mechanism of light nuclei production, the experimental data on particle yields is analyzed in three different setups: (a) no clusters -assumes that light nuclei are formed after the chemical freeze-out, thus light nuclei are to be omitted from the thermal model particle list. Those protons that later-on bind into the light nuclei are counted as 'free' protons at the chemical freeze-out. This is the scenario that had been considered in Ref. [1].
(b) clusters included -light nuclei are formed at the chemical freeze-out. Only stable light nuclei are included in the thermal model particle list.
(c) clusters and decays of unstable nuclei are included -the thermal model additionally includes the feeddown contributions from the decays of unstable A = 4, and A = 5 nuclei to the final yields of protons, deuterons, tritons, 3 He, and 4 He at the chemical freeze-out, as discussed in [31].
The resulting thermodynamic parameters obtained within these three scenarios are presented in Fig. 1 and Table II and compared with the analysis of the data from Ref. [1]. The corresponding data/model ratios for the fitted yields are presented in Fig. 2.
The first scenario, where light nuclei are not included in the particle list, shows a similar fit as the one in Ref. [1]. However, a closer look at the resulting χ 2 profiles in Fig. 1 reveals that this scenario actually has two degenerate solutions. The data can be described by two distinct sets of thermodynamic parameters with similarly good accuracy. The χ 2 profile as a function of the freezeout temperature has two minima located at T ≈ 50 MeV and T ≈ 70 MeV These minima are referred as "low temperature" and "high temperature" minimum, respectively. Both these two minima describe the data well, characterized by χ 2 /N df < 1. Only the T ≈ 50 MeV minimum was discussed in [1] 1 .
To elaborate on the differences between the two minima a detailed look at the role of resonance feeddown is taken. The bottom panel of Fig. 1 (b) shows the fraction of charged pions coming from resonance decays as a function of the temperature for the "no clusters" scenario. It is clear that the feeddown fraction is significantly larger at the high-temperature minimum. This implies that the two minima give similar final yields of measured particles but correspond to significantly different chemical composition of the primordial hadrons. At the high temperature more hadrons are present in the form of excited states while the low temperature system corresponds to essentially a gas of ground state hadrons. 2 Thus, determining the feeddown fractions is one possible way to distinguish the two minima. This can potentially be achieved via a statistical model analysis of yields of short-lived resonances like ρ 0 or K * . Such an analysis may require incorporating partial chemical equilibrium into the HRG model, as discussed in [33]. The yield of the unstable φ meson is already measured by the HADES collaboration [16]. This data may be used in the mentioned partial chemical equilibrium approach, however a careful treatment of the mesons hidden strangeness should be carried out. A strangeness-canonical approach was found to describe well the φ meson yield in lighter systems [34], there a fit to √ s N N = 2.61 GeV Ar+KCl data resulted in the freeze-out temperature of T = 70 ± 3 MeV. Recently it was also pointed out that a thermal model analysis of the unstable ∆-baryon yield in Au+Au collision at the SIS18 energies provides T ≈ 70 MeV which is close to the high temperature minimum [35].
Another difference between the two solutions lies in their thermodynamical properties: the high temperature minimum corresponds to the freeze-out baryon density of n highT B ≈ 0.22 fm −3 , while for the low temperature minimum it is much lower at n lowT B ≈ 0.01 fm −3 . The high temperature description of the fireball created at SIS18 energies may also improve the Siemens-Rasmussen description of the rapidity distribution presented in [1], where the calculated distributions were found to be too narrow as compared to experimentally measured dN/dy.
Both minima describe every measured hadron yield remarkably well, thus, the presently available data alone does not allow to favor one set of parameters over another. Additional information is required. In Sec. IV this question is discussed by utilizing transport model simulations of Au-Au collisions at SIS18 energies.
Next, the possibility that the yields of light nuclei are fixed at the chemical freeze-out is investigated. For this scenario the measured yields of light nuclei are included in the fit and, therefore, the yields of bound protons are excluded from the free protons yield at the chemical freeze-out. Here the fits reveal only a single χ 2 minimum. However, the fit quality worsens significantly, giving χ 2 /ndf > ∼ 20. As was discovered already in 80's [22,36,37] the unstable nuclei play an important role in temperature extraction from the heavy ion collisions. Recently it was pointed out that at lower collision energies the feeddown contributions of unstable A = 4 and A = 5 nuclei to the final yields of protons, deuterons, tritons, 3 He, and 4 He is significant and, for the latter three, may account for as much as 70% of the final yield [31]. The same mechanism as presented in the latter paper is employed here, extending the particle list by 25 unstable nuclei that decay after freeze-out and feed into the final yields of stable light nuclei and free nucleons. The contributions of these light nuclei improve the quality of the fit, the improvement mainly attributed to a better description of the 3 He and 3 H yields. However, the overall description is still poor, with χ 2 /N df ≈ 12.
These results suggest that either the thermal model de-scription of light nuclei at the chemical freeze-out needs to be significantly improved or that thermal nuclei production is not the correct mechanism in the few GeV collision energy regime. Another observation from the fits including the light nuclei is that the description of strange particles is significantly worse than in the no clusters scenario. In general, the fits in all the studied scenarios require values of γ S significantly smaller than unity, meaning that the strangeness is significantly undersaturated relative to chemical equilibrium in the grand canonical ensemble. More involved modeling of strangeness may thus be warranted. In particular, the effect of canonical strangeness suppression may be particularly relevant in this energy range [38].

A. UrQMD model setup
In order to understand the characteristics and dependencies of the fit results from microscopic point of view, a transport model will be used to simulate heavy ion collision data. This allows for the study of the apparent origin and systematic dependencies of the two minima.
Here the heavy ion collisions are simulated with the microscopic transport model UrQMD [39,40], where the yields of all stable hadrons can be extracted at any time step during the evolution. To simulate the most central Au+Au collisions at √ s N N = 2.4 GeV the UrQMD model is used with a restriction b < 4.7 fm on the impact parameter, corresponding to the 10% most central events.
All hadron yields are evaluated assuming the full detector acceptance, and include feeddown from resonance decays. The spectator nucleons are not included in the final yields. In this analysis, no mechanism for light nuclei production is incorporated, hence the UrQMD analysis here corresponds to the "no clusters" setup from the previous section. The UrQMD model used here is based on version 3.4, which is extended here to include an up-to-date set of resonance branching ratios. This extension is essential for a proper description of the sub-threshold strange particle production [41,42]. Nuclear interactions play an important role at the SIS18 energy of the HADES experiment. Therefore, our simulations incorporate the density-dependent nuclear Skyrme potentials. This description gives good results for flow-observables [43] as well as for the space-time evolution of the density [44] at the energies available at SIS18. Figure 3 shows that UrQMD provides a decent description of particle yields measured by the HADES collaboration without invoking any additional free parameters.  denotes thermal fit to all stable particle yields, "UrQMD measured yields" denotes fit to the UrQMD predictions for hadron species that are detected by the experiment. The UrQMD data correspond to particle yields calculated in 4π 10% most central Au+Au events at √ s N N = 2.4 GeV. The UrQMD data correspond to particle yields presented in Fig. 3.

B. Thermal fitting of the UrQMD yields
The UrQMD hadron yields listed in Fig. 3 can be used in a thermal model analysis, in particular to investigate the possibility of the double minimum structure discussed in the previous section. Since the errors of the simulated yields are purely statistical and can be made arbitrarily small, for the thermal fitting procedure a 10% relative systematic error is added, to make it comparable in magnitude to the experimental uncertainties. In this way the thermal fitting procedure for the UrQMD yields is con-sistent with the procedure applied to experimental data in the previous section, thus direct comparisons can be made.
Two possibilities for the set of hadrons included in the fitting procedure are taken: (a) The fit is done only to the hadron yields which are present in the HADES collaboration data (π +,− , K +− , p, Λ).
(b) Yields of all long-lived hadrons from UrQMD are included in the fit (see Fig. 3 for the complete list).
Here are considered separately two options regarding whether the yield of η meson should be included or not.
The χ 2 -profiles of these fits are shown in Fig. 4, where they are also compared with the χ 2 -profile of the fit performed to the HADES data. The fit to UrQMD data that includes only the experimentally measured set of yields shows the two minima structure. This is similar to the finding of the previous section using the HADES data, even the locations of both minima are similar. However, when all the long-lived hadrons, including the η-meson, are included in the fit to the UrQMD predictions, the low temperature minimum becomes rather shallow and then the high temperature fit provides a significantly better description. This improved description may be considered as a direct result of the increased number of independent inputs in the fit, which constrains the fit parameters better. The fit quality decreases slightly, indicating tension of the thermal model with UrQMD. This may be caused by an incomplete chemical equilibrium as well as canonical suppression that can be present in UrQMD.
It is found that the η meson plays a particular role in the analysis. On the one hand, the η meson is a comparatively long-lived particle, such that it decays outside the fireball and thus treated as stable particle in UrQMD.
On the other hand, its lifetime is sufficiently short such that it decays before reaching the detector. One can thus argue whether its yield should be included in the thermal fit or not. If the η meson is omitted from the fit, the χ 2 profile changes such that the low temperature minimum provides a better description than the shallow, high temperature minimum. This suggests that the inclusion of only selected additional particle multiplicities will not necessarily allow for a better distinction between the two χ 2 minima. If the best fit depends strongly on the choice of included hadron multiplicities, these results point to the possibility that the system may have not undergone a universal chemical freeze-out.

C. Time evolution of the solutions
To get a better physical understanding of the meaning of the two minima obtained in the thermal fitting the time evolution of the particle yields in UrQMD is studied following a method suggested in Ref. [45]. To calculate the UrQMD yields corresponding to a given time moment t, the transport evolution is stopped when that time moment is reached and all unstable hadron species are forced to decay. The calculated hadron yields are then related to the respective time t. These calculated yields are then fitted with the thermal model and the extracted thermodynamic parameters are assumed to describe the colliding system at the time t. Here only the experimentally measured yields are analyzed. The time dependence is studied in a range t = 1-20 fm/c.
The results of the thermal fits to the time dependent yields are depicted in Fig. 5 by red (high temperature) and blue (low temperature) bands. The width of the bands is determined by the width of the minima in the χ 2 /ndf surface. The double minimum structure in the fit is found throughout the whole evolution. However, the behaviors of the thermal parameters of each minimum are opposite. The high temperature minimum decreases in temperature and chemical potential (the red band in figure 5), this solution coincides with a "classical" picture of the fireball evolution in a heavy ion collision: temperature and density (chemical potential) decrease with time as a result of the fireball cooling during the expansion. The low temperature minimum on the other hand shows as a function of time a strong decrease of the total entropy per baryon, ∆S/A(t) < 0. At first sight, this seems to violate the second law of thermodynamics which would indeed be a major discovery. However, be reminded that this results from global instantaneous chemical equilibrium fits to a microscopic spatio-temporal non-equilibrium model, which actually respects ∆S/A(t) > 0 throughout.
The results for the time dependence of the temperature are compared with one obtained in a different way, namely via a coarse-graining approach as in Ref. [44] (the green line in Fig. 5). Only the temperature corresponding to the high-temperature minimum is in qualitative agreement with the coarse-grained approach, indicating that the high-temperature is the only physical solution. Note that the agreement is not fully quantitative. In particular the coarse grained temperature, which is extracted from the local densities by thermodynamics relations rather than multiplicities, appears always lower than the chemical temperature. This may be another indication that the system rapidly falls out of chemical equilibrium or may not reach it fully.
Even more important, the entropy per baryon S/A, which can be readily calculated by the HRG, has very different behavior for the two minima as shown in Fig. 5 (d). The entropy along the high temperature minimum has a slight increase, its values of S/A 3 − 5 being close to the values expected for this collision energy [46]. On the other hand, the entropy per baryon of the low temperature minimum behaves abnormally, showing a decrease with time. The values of S/A ∼ 7 − 11 at low temperature appear to be too high for this collision energy. This pathological behavior of the entropy indicates that the low temperature minimum is likely unphysical.
Finally, the corresponding trajectories of the two solutions in the µ B -T plane are shown in Fig. 6 where the respective time is indicated by color. The high temperature represents a trajectory similar to an isentropic expansion [46] for that specific beam energy, while the low temperature evolution resembles a compression with simultaneous heating and loss of entropy.

V. SUMMARY
The hadron yields measured by the HADES collaboration at SIS18 energies [1] reveal the drawback of the thermal model approach which shows ambiguous solutions at this collision energy. This ambiguity is reflected in the existence of two degenerate solutions of the thermal fit to the hadron yields measured in most central Au+Au collisions at √ s N N = 2.4 GeV by the HADES collaboration. While the final hadron multiplicities are almost identical for both statistical model descriptions, the chemical composition and thermodynamic properties show clear differences. The high temperature solution shows a significant contribution of resonance decays to the pion yield, while the low temperature minimum corresponds to a system which consists mainly from ground state hadrons. The high temperature solution is consis-tent with estimates of the ∆ contribution [35] at SIS18 energies, indicating that the higher chemical freeze-out temperature can be confirmed experimentally. The role of light nuclei production at the chemical freeze-out was also studied. The inclusion of these nuclear clusters in the fit significantly worsens the quality of the fit with χ 2 /ndf > 10. In this case only one χ 2 minimum is present with T ≈ 70 MeV, consistent with the high-temperature minimum obtained without the inclusion of light nuclei. The fit improves if the additional feeddown from decays of unstable nuclei is included although the overall data description remains unsatisfactory.
It was found that the particle multiplicities calculated by the UrQMD transport model exhibit the same feature, i.e. can be described by the thermal model in a twofold way, with similar sets of thermodynamical parameters. The detailed study of the time evolution of same multiplicities in the UrQMD model suggests that only one of these solutions, namely the high temperature solution, shows physically reasonable behavior. The low temperature solution, on the other hand, behaves unphysically, with entropy that decreases with time, accompanied by increasing temperature and density. It is concluded that only the temperature of T ≈ 70 MeV may be considered as chemical freeze-out temperature for most central √ s N N = 2.4 GeV Au+Au collisions. This high temperature solution results from a hadrochemical freeze-out that is similar in outcome to the picture of high energy collisions at SPS, RHIC, and LHC. The primordial inelastic nucleon collisions swiftly create a near chemical equilibrium early hot and dense fireball state, which subsequently expands isentropically.
Additionally, it was shown that the inclusion of further, yet unmeasured, multiplicities of stable hadrons would not allow discerning these degenerate states. However, inclusion of unstable hadronic species like K * , ρ 0 , φ, ∆ [33,35] in the fit would seem helpful to resolve the ambiguities. The yield of the unstable φ-meson is already available from the HADES collaboration [16]. This data may be used in partial chemical equilibrium approach [33] to study in more detail the freeze-out mechanism in Au+Au collisions at SIS18 energies, this, however, requires a careful treatment of the quark content of the φ-meson.