Influence of hadronic resonances on the chemical freeze-out in heavy-ion collisions

Detailed knowledge of the hadronic spectrum is still an open question, which has phenomenological consequences on the study of heavy-ion collisions. A previous lattice QCD study concluded that additional strange resonances are missing in the currently tabulated lists provided by the Particle Data Group (PDG). That study identified the list labeled PDG2016+ as the ideal spectrum to be used as an input in thermal-model-based analyses. In this work, we determine the effect of additional resonances on the freeze-out parameters of systems created in heavy-ion collisions. These parameters are obtained from thermal fits of particle yields and net-particle fluctuations. For a complete picture, we compare several hadron lists including both experimentally discovered and theoretically predicted states. We find that the inclusion of additional resonances mildly influences the extracted parameters -- with a general trend of progressively lowering the temperature -- but is not sufficient to close the gap in temperature between light and strange hadrons previously observed in the literature.


I. INTRODUCTION
Determining the number of hadronic resonances and their corresponding properties has been a fundamental question in nuclear physics for decades. Experimentally, new resonances (and their decay channels) are measured through spectroscopy and an up-to-date compilation of those searches can be found in the Particle Data Group booklet [1]. Each new particle is assigned a star on a scale of * to **** which indicates the experimental confidence of the measurement, where * states are the most uncertain and **** states are the most established. Theoretical calculations, such as Quark Models [2,3], predict many more resonances beyond those experimentally measured. Furthermore, there have been general theoretical arguments, originally made by Hagedorn, that near a limiting temperature (now understood as the phase transition of the hadron gas into the Quark Gluon Plasma [4]) an exponentially increasing mass spectrum [5] is expected, i.e. heavier and heavier resonances can form that decay almost immediately into their daughter particles. Using these so-called "Hagedorn states" and comparing thermodynamic quantities to lattice Quan- * Corresponding author: jmstafford@uh.edu tum Chromodynamics (QCD) calculations, initial evidence was found that the older versions of the PDG were insufficient to describe the lattice results [6][7][8]. With the particle production results published by the AL-ICE Collaboration (Pb − Pb 2.76 TeV [9,10] and preliminary Pb − Pb 5.02 TeV [11]) a new issue arose, the socalled "proton anomaly" wherein thermal model fits generally over-predict the proton yields while simultaneously under-predicting the yields of strange baryons. Fluctuation analyses show the same trend [12][13][14]. A number of possible solutions were suggested such as a flavor hierarchy of freeze-out temperatures between light and strange hadrons [15], further missing Hagedorn states [16][17][18], final state reactions [19,20], or corrections to the ideal HRG model [21,22].
In order to settle much of these debates on the number of hadron resonances, lattice calculations have played a vital role in recent years. Initial results from Ref. [16] found evidence of missing strange resonances from first principle calculations. In a more recent work, using partial contributions to the QCD pressure of particles grouped according to their quantum numbers, the exact nature of the missing resonances was determined [23]. Evidence for the need of additional states was seen in certain sectors and even with their inclusion, the contribution of strange mesons still underestimates lattice arXiv:2002.12395v1 [hep-ph] 27 Feb 2020 QCD results. However, in Ref. [23] it was shown that the Quark Model predicts too many states -in particular in the |S| = 2 sector -and thus overestimates most of the partial pressures. The list which correctly reproduces the largest number of observables is the PDG2016+ list (see definition below). For completeness, in this manuscript we will also consider Quark Model states, in order to investigate the effect of including the largest possible number of resonances in the analyses.
It is worth mentioning that missing hadronic resonances have been found to affect a number of areas in heavy-ion collisions such as the transport coefficients η/s and ζ/s [8,24,25], and may somewhat affect collective flow [26,27] and mildly improve the χ 2 of thermal fits [28]. In Ref. [27] we incorporated these missing states from the PDG into a hydrodynamic model and found that they are able to improve the fits to particle spectra and p T . Similar results were found in Ref. [29].
In this follow-up paper, we compare different hadronic lists with an increasing number of states, performing fits of particle yield and fluctuation data from RHIC and the LHC. One of the most pressing questions is whether these new states can close the gap between light and strange freeze-out temperatures, as suggested in Ref. [16]. Initial explorations in this direction can be found in Refs. [30,31], where smaller particle resonance lists were used. We study a scenario with a freeze-out temperature common to all species, as well as one where different freeze-out temperatures result for strange and light species. In the case of fits to yields and their ratios we compare the fit quality from both scenarios.
This manuscript is organized as follows. In Section II we describe the different hadronic lists we utilize in this work. Moreover, we provide the details on how the decay properties of the resonant states were estimated, when not known experimentally. In Section III we briefly describe the HRG model and the data sets employed in the fits we perform. In Sections IV and V we show our results for fits to yields and ratios in the single and double freezeout scenarios, respectively. In Section VI we study the chemical freeze-out from an analysis of net-particle fluctuations from the Beam Energy Scan (BES) at RHIC. We first analyze fluctuations of net-proton and net-charge, then those of net-kaon. Finally, in Section VII we present our conclusions.
The particle lists PDG2016, PDG2016+ and QM utilized in this work, together with the corresponding decay channels can be downloaded from the link provided in Ref. [32].

II. DIFFERENT HADRON LISTS
We consider the following lists of hadronic resonances. PDG2005: PDG from 2005 [33]; PDG2016: PDG **-**** states from 2016 [34]; PDG2016+: PDG *-**** states from 2016 [34]; and QM: PDG2016+ combined with additional Quark Model states from Refs. [2,3]. The com-parison between PDG2005 and PDG2016 is performed to demonstrate the difference between earlier thermal fit models and modern day ones. The PDG2016+ is tested in thermal fits in this paper for the first time. One primary reason why this has not yet been done in the past is that many of the * and ** states do not have adequate decay channel and branching ratio information in order to describe all the particle interactions based only on experimental data. In this paper, we use phenomenology to fill in these gaps. For completeness, we consider a mixture of the PDG2016+ with the addition of Quark Model states that in Ref. [23] appeared to fit some lattice QCD data well. Note that the list we label as QM does not exclusively contain states predicted by the Quark Model. We replaced the Quark Model states that are already known experimentally with the corresponding ones from the PDG2016+ list. In order to discern whether a state predicted by the Quark Model was already listed in the PDG, we considered states sharing the same quantum numbers (baryon number, strangeness, electric charge, isospin and spin degeneracy), with comparable masses (mass difference below ∼ 150 MeV, or within 15% of the mass of the PDG state, whichever is smaller). The issue with the decays is even more pronounced for the QM, where no decay channels are given. In this case we relied entirely on phenomenological approaches that are described in detail below.

A. Decay lists
In the PDG database [34], only the relatively few established states are provided with an exhaustive set of decays, whereas most states are only given branching ratios which do not sum up to 100%. In this case, additional decays are listed as "dominant", "seen" or "not seen".
When considering the effect of resonance decays, we include all the strong decays listed by the PDG with a branching ratio of at least 1%, but we do not include weak decays. When the listed branching ratios do not sum up to 100%, we assign the remaining percentage to an electromagnetic decay of the kind N 2 → N 1 + γwhere N 2 and N 1 are hadrons with the same quantum numbers, N 1 being the next state in descending mass order with compatible parity for such a decay. We decay resonances that lack quantitative information predominantly radiatively (as explained above), in order to populate the next state with known branching ratios.
The situation is different for the hadronic resonant states predicted by relativistic quark models, since these are not supplied with any information on their decay properties. In order to take into account the effect of these additional states on observables, we need to estimate their branching ratios.
Along the lines of the work done for the PDG states, we infer the branching ratios of the additional QM states using the following procedure. Every state is assigned a predominant radiative decay, analogously to PDG states with incomplete information. The remaining percentage is reserved for hadronic decays, for which we use linear extrapolation from known branching ratios of PDG particles with the same quantum numbers B, S, I and Q (baryon number, strangeness, isospin and electric charge).
The linear extrapolation is performed as follows. States from the PDG2016+ list are divided into "families" with the same quantum numbers, and all strong decay modes present in the family are grouped together. For each channel appearing in the family, a linear interpolation of its branching ratio as a function of the particle mass is performed. Then, QM states of the same family are assigned branching ratios by sampling the linear mass dependence constructed as explained above. At this point, all negative BR values are discarded, as well as all decays that violate mass conservation. Finally, the sum of branching ratios for strong decays is normalized and put together with the electromagnetic ones.
We emphasize here that for the QM states we have the least amount of information and, therefore, the largest amount of uncertainties.

III. THERMAL FITS AND THE HADRON RESONANCE GAS MODEL
The study of chemical freeze-out we present here makes use of the Hadron Resonance Gas (HRG) model, which describes the hadronic phase below the transition temperature T c as a system of non-interacting particles [5,35,36]. The HRG model has a wide-spread applicability in heavy-ion studies in reproducing thermal particle abundances and lately also in providing results on fluctuations of conserved charges (B, Q, S) [12,13,[37][38][39]. Recently, these observables have been compared to the measured moments of net-particle distributions, and provided freeze-out temperatures which are compatible with the ones obtained by comparing the experimental data to lattice QCD calculations [37,[40][41][42][43][44].
Historically, the HRG model has been widely employed to compare data on particle production for energies ranging from the AGS to the LHC [45][46][47][48][49][50][51]. Produced particle yields N i are obtained by adding the contribution from resonances to the primordial thermal yield, given by V n i : In the above, n i R is the average number of particles of type i resulting from a decay of resonance R, n i and n R are thermal densities calculated through the statistical model, and V is the system volume. The decay of resonance R into stable particles such as p, π, K, Λ, Ξ and Ω is taken into account by introducing an effective chemical potential µ R = i µ i n i R . Here µ i = B i µ B +Q i µ Q +S i µ S is the chemical potential for particle i carrying specific baryonic, electric charge and strangeness quantum numbers. Conditions on the net-strangeness and net-charge density are imposed, to match the heavyion collision situation: These allow one to constrain the three chemical potentials. In this way, yields and ratios calculated within the HRG model only depend on the thermal parameters (T, µ B ) (and V in the case of yields). Thermal properties at the chemical freeze-out have been studied using yields and ratios from STAR data in Au-Au collisions at √ s NN = 200, 39, 27, 19.6, 11.5, 7.7 GeV [52][53][54] and from ALICE data in Pb-Pb collisions at √ s NN = 2.76 TeV and 5.02 TeV [9,11,[55][56][57].
In this manuscript, we focus on the STAR data at √ s NN = 200 GeV and 0 − 5% centrality, and the AL-ICE data at √ s NN = 5.02 TeV and 0 − 10% centrality.
We perform thermal fits of the particle yields and ratios, using published data on ratios, if available, for STAR [52] and for ALICE [9]. In order to build the remaining ratios, published data on yields from both collaborations have been used with a proper propagation of the errors in the final result. We evaluate the yields and ratios for each hadronic list and extract the thermal parameters (T, µ B , V ) by using the thermal fit package FIST [58]. The package allows users to choose their own particle lists, as well as data sets, in the fit.

IV. SINGLE FREEZE-OUT SCENARIO
Initially, we consider a common freeze-out temperature for strange and light hadrons.
We fit the measured yields for the following particles: π + , π − , K + , K − , p,p, Λ,Λ, φ, Ξ,Ξ, Ω,Ω at the LHC, while at RHIC the separate Ω andΩ yields are replaced by the sum Ω +Ω. When we take the ratios, we divide the light particle yields by the yield of either π + or π − and the strange particle yields by the yield of either K + or K − . The results of the thermal fits for both yields and ratios while varying the particle resonance list are shown in Table I for LHC data at √ s N N = 5.02 TeV and in Table II for STAR data at √ s N N = 200 GeV. At the LHC we hold µ B = 1 MeV fixed to avoid the possibility of negative chemical potentials and remain consistent with previous analyses [59]. From Tables I and II a few trends begin to emerge: • for both yields and ratios, more hadronic states generally decrease the chemical freeze-out temperature; • the chemical freeze-out temperatures from yields and ratios approximately agree; • generally, due to the decrease of the chemical freeze-out temperature with the increase of the number of hadronic states, the volume increases accordingly; • µ B is relatively stable when changing the particle resonance list; • RHIC appears to have somewhat higher chemical freeze-out temperatures compared to the LHC; • the inclusion of resonant states beyond the PDG2016+ list improves the fit quality at the LHC for both yields and ratios, while at RHIC the PDG2016+ fits best for the yields and the QM and PDG2016+ are equivalent for the ratios.
The second to last point is puzzling because one expects the freeze-out curves to remain roughly constant at low baryon chemical potential, according to lattice QCD studies of the phase transition at finite µ B [60,61].
The values of temperature and chemical potential listed in Tables I and II are shown in Fig. 1 and Fig. 2, respectively.
In Fig. 1, the decreasing trend of the chemical freezeout temperature with the increasing number of hadronic states in the list is clearly visible. As mentioned previously, slightly higher temperatures are found at RHIC. On the other hand, in Fig. 2 we can see that the RHIC freeze-out chemical potential barely changes with different resonance lists, and that the results from yields and ratios agree well.

V. TWO FREEZE-OUT SCENARIO
We now explore the possibility of a flavor hierarchy between light and strange particles, using the same particle resonance lists as in the previous section. There,   IV. Temperatures and baryon chemical potentials in the double chemical freeze-out scenario obtained from fits to ratios or total yields from STAR data [52,53] for 0 − 5% centrality in AuAu collisions at √ sNN = 200 GeV, using different PDG lists. For the case of total yields, the volume is shown as well. The last two columns show the chi square per number of degrees of freedom in the fits. we saw that the Quark Model list produces the best fit to the experimental data at ALICE, while at RHIC the PDG2016+ provided the best fit for the yields and for the ratios the QM and PDG2016+ list were nearly equivalent. Nevertheless, we stress here again that the QM results are shown here just for illustration, as in Ref. [23] it was found that the QM contains too many states and fails to reproduce several thermodynamic quantities from lattice QCD. For this analysis, we allow for two separate chemical freeze-out temperatures as well as two separate freeze-out baryon chemical potentials for light and strange particles.
For the light hadron fit we use π + , π − , p, p, K + and K − , while for the strange hadrons fit we use all species except pions and (anti)protons. We construct the ratios of strange hadrons dividing by the K ± yields, to avoid mixing of light and strange degrees of freedom. We do this also because, if light particles freeze-out at a lower temperature than the one for strange particles, the respective volumes would be different and thus would not cancel out when taking ratios. In Tables III and IV the extracted chemical freeze-out temperatures, chemical potentials and volumes (for yields only, as before) are shown for both the light and strange hadrons at the LHC and RHIC, respectively. As in the case of a single freeze-out scenario, there is very good agreement between fits to yields and ratios. Moreover, the freeze-out temperature for strange hadrons is systematically higher than that of light hadrons.
An additional effect of the two temperatures is that the light temperature is significantly lower than in the single freeze-out temperature scenario. This means that the fit in that case is driven by the strange states, as it was discussed in [38]. We also notice that the inclusion of more resonances affects the strange freeze-out temperature more than the light one. This is not unexpected, as most of the additional states carry strangeness.
A summary of our results from Tables III and IV is shown in Fig. 3, for yields and ratios from both ALICE and STAR, and with the four different hadron lists we consider. In this pictorial representation, the separation in temperature between strange and light freeze-out appears very clearly, especially from fits to LHC data. In this case, all lists lead to a very pronounced separation. In the case of STAR data, the separation is especially evident when utilizing the PDG2016 and PDG2016+ lists, although still present even for fits with the QM list. There seems to be a general trend in these last three hadron lists, where the presence of more states leads to a smaller separation of the strange and light freeze-out temperatures. The PDG2005 list is excluded from this trend: it contains a significantly smaller number of resonances, which is not realistic in modern calculations. On the other hand, we also know that the strange states contained in the QM list are too numerous, as noted earlier and shown in Ref. [23].
In Appendix A (Figs. 5-8) a comparison of the yields and ratios to experimental data is shown, both in the single (red lines) and double (blue dotted lines) freezeout scenarios.

VI. FREEZE-OUT FROM FLUCTUATIONS
In this Section we complete our analysis by studying net-particle fluctuations from the BES at RHIC. We employ the state-of-the-art list PDG2016+ to perform the analysis of net-proton and net-charge fluctuations from Ref. [37], as well as the analysis of net-kaon fluctuations from Ref. [12]. Both analyses were originally carried out with an older PDG list, which we here indicate as PDG2012.
In order to determine both the temperature and chemical potential at the chemical freeze-out, in general two experimental quantities are required. However, due to the large experimental errors on higher order fluctuations (i.e., the ratio χ 3 /χ 2 ) as noted in [62,63], in Ref. [12] we obtained the freeze-out parameters as follows. Starting from the net-proton and net-charge freeze-out parameters obtained in Ref. [37], we followed the isentropic trajectories determined in Ref. [64] via a Taylor-expanded equation of state from lattice QCD. These were constructed by first determining the entropy-per-baryon ratio at the freeze-out point for each collision energy from Ref. [37], and then following the path that conserves S/N B 1 . The ratio χ K 1 /χ K 2 was then calculated along these trajectories , and compared to the experimental 1 Isentropic lines are the trajectories followed by the system in the case where no dissipation is present. This is a good approximation due to the extremely low viscosity of the medium created in results. For each collision energy, the overlap region will then correspond to the freeze-out points for net-kaons, which are shown in Fig. 4 in gray, while the red points correspond to the net-proton and net-charge freeze-out points from Ref. [37]. In order to repeat the analysis with the PDG2016+ list, we first calculated the freeze-out parameters for light hadrons via a combined fit of χ p 1 /χ p 2 and χ Q 1 /χ Q 2 with the new list. We then used these light freeze-out parameters to determine the new isentropic trajectories in the QCD phase diagram. As was done in Ref. [64], we utilized a Taylor-expanded equation of state from lattice QCD. With the new isentropic trajectories, we could then calculate the net-kaon fluctuations along them and determine the freeze-out parameters.
We note that, both for the net-proton and net-charge, as well as for the net-kaon analysis, we included the effect of resonance decays, and imposed the same kinematic cuts employed in the experiments [65][66][67]. Moreover, for the light hadron analysis we included the effect of isospin randomization [68,69], which has a large impact on netproton fluctuations.
The results for the light and kaon freeze-out with the PDG2016+ list are shown in Fig. 4 with black points and blue shapes, respectively. We can see that, as in the analysis with the PDG2012 list from Refs. [12,37], there is a separation of the net-kaon from the light particle freeze-out temperature at the higher collision energies (lower µ B ) that decreases as the collision energy decreases. With the inclusion of more resonances in the HRG model, the freeze-out temperature for the net-kaons is decreased, while the light hadron freeze-out is mostly unchanged. However, the gap between the two freezeout temperatures is not closed, and the resulting temperatures are compatible with the ones obtained in the previous Section with the analysis of fits in the double freeze-out scenario.

VII. CONCLUSIONS
In this paper we explored the effects that including additional, not-well-established hadronic resonances in the PDG lists, as well as a combination of the latter with predicted (but not yet measured) Quark Model states, has on the extracted freeze-out parameters in heavy-ion collisions. In Ref. [16] it was suggested that the addition of missing strange states could possibly explain the tension between proton and multi-strange hadrons production at the LHC. An alternative picture was proposed in Ref. [15] wherein strange hadrons freeze-out at a higher temperature than the light ones, which is due to the higher hadronization temperature of strange particles heavy-ion collisions, and more so if they are utilized in a small portion of the evolution, in the proximity of the transition, as we do here.   compared to light particles. We explore this alternative picture here.
For the first time in this paper, we systematically ex-plore the effect of enlarged hadronic resonance spectra -along with their decay properties -on thermal fits, and find that this tension is not resolved. Especially at the LHC, the separation between strange and light freeze-out is extremely pronounced. At RHIC the separation is smaller, and is almost resolved when states from the Quark Model are included, although it is important to remember that certainly too many strange states are predicted by these calculations, as was shown in Ref. [23]. From thermal fits based on the most realistic particle list PDG2016+, we estimate the light chemical freeze-out temperature to be T L ∼ 138 − 144 MeV and the strange chemical freeze-out temperature to be T S ∼ 162 − 167 MeV at the LHC, and T L ∼ 148 − 158 MeV and T S ∼ 160 − 165 MeV at RHIC. Generally, the light chemical freeze-out temperature appears to be consistent with the net-p and net-Q results (within error bars), while the strange temperature is consistent with the one obtained from net-K fluctuations.
Finally, we note that in other works [17,18,20] it was suggested the a dynamical scenario could explain the proton to pion puzzle at the LHC. The inclusion of the additional states in such a scenario is left for a future work. Generally, we find hints of two separate chemical freeze-out temperatures, but smaller error bars in the experimental data would be needed before a definitive answer can be given.