Systematics of chemical freeze-out parameters in heavy-ion collision experiments

We discuss systematic uncertainties in the chemical freeze-out parameters from the $\chi^2$ analysis of hadron multiplicity ratios in the heavy-ion collision experiments. The systematics due to the choice of specific hadron ratios are found to lie within the experimental uncertainties. The variations obtained by removing the usual constraints on the conserved charges show similar behavior. The net charge to net baryon ratios in such unconstrained systems are commensurate with the expected value obtained from the protons and neutrons of the colliding nuclei up to the center of mass energies $\sim 40$ GeV. Beyond that the uncertainties in this ratio gradually increases, possibly indicating the reduction in baryon stopping.


INTRODUCTION
Experimental data from the relativistic heavy-ion collisions indicate that a considerable fraction of energies of the colliding nuclei are deposited in the region of interaction to form a hot and dense fireball. The resulting system may thermalize rapidly to form the exotic state of thermally equilibrated quark gluon matter [1][2][3]. Pressure gradients from the central to the peripheral region of the system would lead to expansion and cooling. Thereafter the system would most likely make a transition from the partonic phase to the hadronic phase in thermal and chemical equilibrium [4,5]. With further expansion, and increase in inter-particle separation, the chemical composition of the hadrons would freeze out. The characteristics of the freeze-out surface depends on the nature of interaction or scattering [6].
Several statistical models of hadronic system has been formulated in last few years, based on the seminal work by Fermi [7], followed by Pomeranchuk [8] and Landau [9,10] in subsequent years. Hagedorn recognized that the existence of numerous hadronic resonances is a characterization of hadron's strong interactions [4]. These simple statistical approaches are successful with surprising accuracy in explaining the yield of experimentally observed hadrons produced in collisions of elementary particles and as well as heavy-ions [6,. The underlying consideration in these analyses are that, even when the hadronic interactions cease at some point of evolution of the fireball, the hadrons populate the phase space according to statistical distribution [21].
Assuming the chemical equilibration of the full hadronic spectra, the thermodynamic freeze-out parameters are obtained from the statistical model calculations by performing a χ 2 fit with the available experimental multiplicity data [25-27, 30, 47-53]. One of the most important issues has been the discussions of resulting enhancement of multi-strange hadrons [54], as well as the saturation of strangeness [2]. However possible scenarios of departure of strangeness chemical equilibrium introduced through an additional parameter − the strangeness suppression factor γ s , have been discussed [3,20,24,25,27,52,55]. More recently the possibility of mass dependent or flavor dependent freezeout pictures have emerged (see e.g. [48,53]). The studies in this direction has reached the state-of-theart with several sophisticated codes like THERMUS [56], SHARE [57], THERMINATOR [58] publicly available for studying the freeze-out characteristics.
Here we revisit the analysis of the chemical freeze-out parameters from the χ 2 fit of experimental multiplicity data ratios using the Hadron Resonance Gas (HRG) model for a different purpose. It has been generally found that equilibration parameters are biased to the chosen particle ratios. For example in Ref. [25] a particular set of independent hadron ratios are used with minimum repetition of any given hadronic species to reduce systematic bias. We shall instead consider a variety of possible sets of hadron yield ratios to obtain the variations of the estimated freeze-out parameters, and the resulting variations in the predicted yield ratios. This would give us a quantitative estimate of the systematic uncertainty in the analysis method due to the choice of the hadron yield ratios.
Another possible source of systematic uncertainty is the choice of constraints to be satisfied by the system. The most popular choice is to consider a fixed ratio of net electric charge to net baryon number depending on the colliding nuclei, along with the setting of net strangeness to be zero. Usually the available hadronic yield data is obtained within a limited central rapidity bin. Though the constraints should definitely be satisfied globally, there is no imperative reasons to enforce these conditions in a given rapidity bin. For example, at very high center of mass energy of collisions ( √ s) one expects baryon stopping in the central rapidity region to be negligible. The system may simply be composed of the secondary particles produced, and would have both the net charge and net baryon number to vanish. It is therefore natural to study the system without considering any constraints at all. Here we shall perform this analysis and check the conditions on the net charges when the constraints are not used.
The paper is organized as follows. A brief description of HRG model is given in section II. In section III − VI we have described our model and its results followed by discussion. Our conclusions are presented in section VII.

II. HADRON RESONANCE GAS MODEL
For strong interactions, chemical equilibrium means the equilibration of the conserved charges baryon number (B), electric charge (Q) and strangeness (S). Thus, the equilibrium thermodynamic parameters are the temperature T and the three chemical potentials µ B , µ Q and µ S , corresponding to the three conserved charges respectively. At chemical freeze-out, the HRG model may provide a good description of the system in thermodynamic equilibrium [14,15,[17][18][19][24][25][26]. The partition function of HRG in the grand canonical ensemble is given as, where sum runs over all the hadrons and resonances. Here thermodynamic potential for i'th species is given as, where the upper sign is for fermions and lower is for bosons, and V is the volume of the system. Here g i , E i and m i are respectively the degeneracy factor, energy and mass of i th hadron. The chemical potential of the i th hadron may be expressed as, µ i = B i µ B + S i µ S + Q i µ Q , with B i , S i and Q i denoting its baryon number, strangeness and electric charge.

III. APPLICATION TO FREEZE-OUT
From the partition function, the number density n i is obtained as, (3) One may relate this number density for the i'th detected hadron to the corresponding rapidity density as [27], where the subscript Det denotes the detected hadrons. Here, where the summation is over the unstable resonances j that decay to the i th hadron. To remove systematics due to the volume factor the ratios of yields are usually employed. Thus, the expected equations to be satisfied are, where, R exp α is the ratio of hadron yields obtained from experiments and R therm α is the ratio of the number densities from the model calculations. However individual hadron numbers are not conserved in strong interactions. Rather an equilibrium of the hadrons in terms of the conserved charges B, Q and S are sought. Therefore a χ 2 fitting is performed with as many such ratios as permissible.
It is a common approach to first fix two of the parameters, say, µ Q and µ S using the constraint relations [47], and The value of the ratio of net baryon number to net charge depends on the physical system. For example, in Au + Au collisions, constant = N p /(N p + N n ) ≃ 0.4, with N p and N n denoting the number of protons and neutrons in the colliding nuclei. The problem is now reduced to obtaining the best fit value of the other two parameters by minimizing the χ 2 , which is defined as, Here σ exp α is the experimental uncertainty in the ratio R exp α . This kind of a setup, and its various improved versions in terms of the modifications of the HRG model, has successfully described hadron yields from AGS to LHC energies [14-21, 24-28, 47, 48, 50-53, 59].

IV. SYSTEMATIC UNCERTAINTIES
There are many possible sources of systematic uncertainties that enter in the analysis. The systematic uncertainties in the experimental data are expected to be reduced when considering the ratios of hadron yields. For example, in Ref. [25,55] more particle-antiparticle ratios are used than the cross ratios, as they may achieve some cancellation of the inherent systematic errors of experimental measurements.

A. Systematics due to choice of hadron ratios
However the use of the hadron ratios introduces further systematics into the picture. Therefore different authors consider different criteria for choosing the set of hadron ratios. For example Ref. [25] claimed that whenever a given experimental yield is used several times in the ratios there may be some correlation uncertainties in particle ratios. An estimate of such uncertainties was found to be less than 5 percent Ref. [55]. On the other hand, choosing all possible hadron yield ratios do not serve the purpose either, as there will be relevant correlations between different ratios [25,60]. The alternate possibility is to choose (N − 1) number of statistically independent ratios of N hadron yields. But it has been argued that this procedure may lead to the information loss [27,60].
Here we propose to quantify the systematic uncertainty due to the consideration of different sets of statistically independent hadron yield ratios. With the N experimentally measured hadron yields one can form N C 2 = N (N − 1)/2 yield ratios (interchanging numerator and denominator does not give any new information). Any (N − 1) of these ratios are statistically independent as long as each of the N hadron yields appear at least once. For example, if there are three hadrons π + , π − and p, there are in all 3 C 2 = 3 possible ratios that can be formed. Out of these 3 ratios only (3 − 1) = 2 are independent, because the third ratio can always be written as a combination of the 2 independent ratios. However any 2 of these three ratios can be chosen to be independent, and there are 3 different ways to choose them. For higher number of hadrons, all the sets of (N − 1) ratios may not necessarily be independent. Still the number of independent sets would grow very large with the increase in number of hadron yields. In general, for N given hadron yields there are N (N − 1)/2 possible ratios, from which (N − 1) ratios may be chosen in N (N −1)/2 C N −1 ways.
One can obtain the freeze-out parameters for each of these sets of independent ratios. The variation of the parameters with the choice of these different sets would then be an useful measure of the concerned systematic uncertainty.
For the constrained system, we shall use Eq. 7 and Eq. 8 for fixing µ B and µ S , and fit T and µ Q from the minimizing conditions on χ 2 given as, We shall perform the χ 2 analysis for each of the chosen sets of ratios and obtain the variation of the resulting freeze-out parameters. Finally we shall obtain the corresponding variations of the predicted hadron yield ratios.

B. Relaxing the net charge constraints
The other possible systematics is the use of the constraint equations Eq. 7 and Eq. 8 respectively. The constraints essentially ensure that the overall net baryon to net charge ratio as well as the vanishing of net strangeness is preserved even locally for the small rapidity window of the experimental data. However such restrictions are not governed by any physical laws. In fact for very large √ s, say for LHC energies, baryon stopping is expected to be negligible. In that case both net baryon number and net charge may vanish and their ratio may be indeterminate. Therefore one can try to obtain the freeze-out parameters by varying the constraints. For example in Ref. [55] variations in the systematics were obtained by choosing several conditions like (i) µ Q = 0, (ii) µ Q constrained by net baryon to net charge ratio, and (iii) µ Q constrained to that ratio along with vanishing net strangeness. Similarly there has been studies [61] that discuss possibilities of net strangeness being non-zero in these regions. The most drastic change would be to remove the constraints altogether, which we shall try here. The thermodynamic parameters will be obtained from the four χ 2 minimization equations, Here also we shall study the systematics due to the choice of hadron yield ratios as discussed in the last subsection. Further it would be interesting to study the systematics of the net baryon to net charge ratio and the net strangeness with varying √ s.

V. DATA ANALYSIS
We have used AGS [62][63][64][65][66][67][68][69][70], SPS [71][72][73][74][75][76][77][78][79][80], RHIC [81][82][83][84][85][86][87][88][89][90][91][92][93][94][95][96] and LHC [97][98][99][100] data for our analysis. The data for STAR BES is obtained from Ref [48,55,101]. The data at mid-rapidity for the most central collisions are considered. All hadrons with mass up to 2 GeV are included for the HRG spectrum. The masses and branching ratios used are as given in [56,102]. Experimental yields are available for only a few hadrons at various collision energies. Also the data for all the hadronic yields at one collision energy may not be available in another. For example, the LHC data set does not have theΛ, and we assumed it to be same as that reported for Λ. Similarly, we did not use Ω data for any parametrization, as the individual yields of Ω + and Ω − are not available for most of the √ s. Data for φ (1019.46 MeV) was not used in the fitting as it is already included in the analysis through its strong decay channel to kaon. So the yields of the hadrons used to obtain the freeze-out parameters are, π ± (139.57 MeV), k ± (493.68 MeV), p,p (938.27 MeV), Λ,Λ (1115.68 MeV), Ξ ∓ (1321.71 MeV). For the AGS at √ s = 4.85 GeV, Ξ data was not available. For the lower AGS energies we used π ± , k ± , p and Λ data. For the ten hadron yields in SPS, RHIC, and LHC data, one can construct a set of nine independent ratios as discussed earlier. The ratios in each set will be independent if each of the ten hadron yields appear at least once in the nine ratios. Therefore not all choices of nine ratios will contain independent ratios. There is a huge multiplicity of such suitable sets. We selected a sample of only nine such sets as given in table I, to make an estimate of the systematic uncertainty due to the choice of sets. These sets were chosen such that each of the 10 C 2 = 45 ratios appear in at least one set. However some of the ratios like π − /π + appear in all of them and some others appear less. Even some of the ratios like k + /π + and Ξ + /Ξ − appear only once. We shall extract freeze-out parameters for each of the given sets. The variation of the parameters from these sets will be our estimate for the systematic uncertainties.
The central values would be obtained from the weighted mean. Investigating the effects of choosing few more such sets may be attempted in future, but we note that the actual number of all such independent sets is prohibitively large.
All these considerations will remain same when we study the freeze-out characteristics by removing the constraints on the net charges. Our numerical code solves the equations Eq. [7 − 11] for the case including the constraints, and the equations Eq. [12− 13] for the case without the constraints using the Broyden's method with a convergence criteria of 10 −6 or better.

A. Freeze-out Parameters
The variation of the extracted freeze-out parameters with the center of mass energy √ s is presented in Fig. 1. The error bars indicate the variations of the parameters obtained from the different sets of hadron yield ratios. Results including the constraints on net charges and those for the unconstrained systems are shown separately. We find that the extracted parameters in both the schemes are in reasonable agreement with each other. For comparison, the freeze-out parameters in some of the recent literature are also plotted. In general our fitted parameters are commensurate with those reported in literature with the various forms of the HRG model and also the various choices of the computation technique. For example Ref. [20,24,25,27,52,55]  ered a strangeness suppression factor. For small value of µ Q , it is considered to be zero in Ref. [55]. In another work, the fireball volume was used as a parameter in Ref. [25,26,30,48]. Interestingly the error bars considered in the literature, which are obtained by varying the χ 2 by a small amount, agree well with the systematic variations obtained by us. In Fig. 1a, the freeze-out temperatures T are shown. The temperature rises with increasing √ s and approaches a saturation [4], except at the LHC energy where the temperature is lower [100,103].
The freeze-out chemical potentials µ B , µ Q and µ S , as functions of √ s, are shown in Fig. 1b. The baryon chemical potential is supposed to be larger for lower √ s due to baryon stopping. This is also why a non-zero charge chemical potential arise due to the net isospin of the colliding nuclei. On the other hand a non-zero strangeness chemical potential requires different explanation as there is no net strangeness in the colliding nuclei. The non-zero baryon chemical potential induces production of strange baryons. This in turn requires a strange chemical potential, so that enough strange anti-particles are produced to ensure the vanishing of net strangeness. The appearance of a non-zero strangeness chemical potential is perhaps the strongest signal of the chemical equilibrium of the strongly interacting system during its evolution [104].
One would expect such a picture to hold naturally in the scheme where the net charge constraints have been imposed, ensuring the vanishing of net strangeness. However we find that the strange chemical potential to follow the same quantitative behavior even when the constraints are removed. In fact all the freeze-out parameters agree quite well within the systematics. It appears that the system is so close to chemical equilibrium that the inclusion or exclusion of external constraints become irrelevant. The goodness of the fit can be determined by the reduced χ 2 i.e. the ratio of χ 2 with the number of degrees of freedom (NDF), which is shown in Fig. 2 as a function of √ s. We show the comparison between our two schemes of χ 2 minimization and also plot the reported results available in the literature. We find the reduced χ 2 in the two schemes to be generally in agreement with each other. We reiterate that the variations for each scheme are due to the various choice of independent sets of hadron yield ratios. The widest variations observed are for the SPS data set, and for the lowest AGS data set.
As seen from Fig. 2, our results are in close proximity to the results already reported in the literature. It is a general norm to choose a given set of particle ratios so that the reduced χ 2 lie close to 1. However this is not a necessary condition. The only requirement is that if we obtain the reduced χ 2 for the variety of independent sets of hadron yield ratios, they should have a χ 2 distribution with the mean value close to the number of degrees of freedom. Therefore it should not be considered alarming if the reduced χ 2 value is far from 1. With the few sets of ratios chosen we find the mean reduced χ 2 to be close to 1 for most of the √ s. A detailed study with a larger number of sets will be reported elsewhere.

B. Hadron Yield Ratios
With the fitted freeze-out parameters we now obtain the various hadron yield ratios and compare them with the experimental data. Here we shall discuss only some of the relevant hadron ratios. The variation of π − to π + ratio with collisional energy is shown in the Fig. 3a. This ratio depends only on T and µ Q . The slight over abundance of neutrons to that of protons in colliding nuclei leads to an isospin asymmetry yielding more π − over π + . Thus the experimental π − to π + yield ratio is greater than 1 for lower √ s, eventually approaching 1 at higher collision energies. The extracted µ Q is therefore slightly negative for small √ s and tends to zero for the higher √ s. The ratios obtained from our parametrization for both the schemes follow the same trend throughout the collisional energy range and reproduced the experimental data quite well.
In Fig 3b, the extracted ratios p/π + andp/π − are similarly found to agree well with the experimental data. These ratios are dependent on the interplay between T , µ B and µ Q . At lower √ s, p >p and gradually approaches 1. On the other hand the π − is always almost close to π + . This leads to the converging trend of the two shown ratios at higher √ s. The same trend is also followed by the k + /π + and k − /π − ratios as shown in Fig. 3c. Apart from T and µ Q , this ratio depends on µ S rather than µ B . However the behavior of the k + /π + ratio is reversed for the lower √ s, giving rise to a ′ horn ′ like structure near √ s ≃ 10GeV. This was originally suggested as an observable for strangeness enhancement, a signature of onset of deconfinement and QGP formation [105][106][107][108]. The explanation of this behavior may well be beyond the scope of the HRG model, but the model predictions agree well with the experimental data.
Dependence of Λ/p and Ξ − /p ratios on collisional energy are shown in the Fig. 4a and Fig. 4b, respectively. The agreement for Λ/p with experimental data is reasonable except for a slight under prediction from model analysis. Consideration of possible uncertainties in contribution from weak decays may remove this discrepancy [89]. However, the model predictions for Ξ − /p is found to agree with the experimental data.
Finally in Fig. 5 we plot the variation of some hadron yield ratios whose experimental data were never used in our analysis. Specifically we choose the multi-strange hadrons Ω and φ. Note that though their yields from the experiments are not used they appear implicitly when we consider their thermal abundance and their contribution in the decay chain of the hadrons. Our predictions seems to be in excellent agreement with experimental data and the systematic uncertainties for both the χ 2 schemes are within control.
This brings us back to the discussion of the bias induced by the choice of the the set of yield ratios for obtaining the freeze-out parameters. As discussed earlier that various authors used suitable sets of hadron ratios so that the value of the χ 2 is close to 1. We discussed above that this is not a necessity. Further, our results show that the estimated systematic uncertainties due to the choice of different hadron ratio sets are all within the experimental bounds. Among the sample ratio sets, while the π − /π + ratio has been used in all the nine sets considered, the other four ratios have only been used in one among those sets. Still the predicted values are close to the experimental data. Obviously this would not have been possible if the system is very far from the equilibrium description. However we should note that only a small subset of all possible sets of independent ratios were actually used in our study.
In a recent work [109], we proposed an alternate scheme to the χ 2 analysis for studying the freeze-out characteristics. The basis of that work was to argue that one should use the quantities associated with various conserved charges of strong interactions to obtain the freezeout parameters. But sufficiently close to equilibrium either of these approaches would produce commensurate results. Our present results with the considered systematics are found to be quantitatively in agreement with that alternate prescription.

C. Effects of removing the net charge constraints
We finally discuss the predicted behavior of the net charges for the unconstrained system. In Fig. 6, we show the variation of NetQ/NetB as well as the NetS as a function of √ s. We find that the NetS remains consistent with zero for the full range of collision energy even without any constraint being put in. It seems that for chemical equilibration, the system is essentially guided by the strangeness sector both in terms of the strange chemical potential as well as the net strangeness. At the same time the NetQ/NetB is also consistent with the expected value of 0.4 up to √ s = 39 GeV. Thus even without any constraint the system seems to obey the relations. If the constraints should hold in the central rapidity region when baryon stopping is significant, then possibly beyond √ s ∼ 40 GeV, baryon stopping is no longer favorable.
For higher √ s the central value of NetQ/NetB seems to decrease, along with a larger systematic uncertainty that increases with increasing √ s. The individual NetQ and NetB are unconstrained inside the central rapidity region so long as these charges are conserved globally. So there can be possible uncertainties in their ratios depending on the chosen set of hadron yield ratios. Also for very large collision energies both of these charges should approach π -/π + √ S NN in GeV π -/π + without constraints π -/π + with constraints Experimental π -/π + (a) p/π -√ S NN in GeV p/π + without constraints p/π + with constraints Experimental p/π + p/πwithout constraints p/πwith constraints Experimental p/π - /π + without constraints k + /π + with constraints Experimental k + /π + k -/πwithout constraints k -/πwith constraints Experimental k -/π -(c) FIG. 3. Variation of π − /π + , p/π + ,p/π − , k + /π + and k − /π − with √ s zero. This would imply that the ratio could fluctuate depending on which of the two charges are estimated to vanish with respect to the other charge. The systematic uncertainty is therefore found to increase with √ s and becomes maximum at the LHC energy.

VII. SUMMARY AND CONCLUSION
Most of the existing literature employ the χ 2 analysis as the method of choice for extracting the chemical freeze-out parameters from the available hadron yields. The usual practice is to choose a suitable set of hadron ratios such that the χ 2 is minimized with respect to the freeze-out parameters close to the value of χ 2 /N DF = 1.
To quantify the sensitivity of such a bias, we performed the χ 2 analysis with various sets of hadron yield ratios. We find that the variation of the extracted chemical freeze-out parameters are very much under control, as the resulting systematic uncertainties for the hadron yield ratios lie within close proximity to the experimental bounds.
Moreover in the standard χ 2 scheme there are two common assumptions. The first one is that of a fixed net charge to net baryon number ratio depending on the colliding nuclei, and the other is that of strangeness neutrality. Although these should indeed be true for the whole phase space of a given collision event, they may not necessarily hold for the central rapidity bin. Here we studied a parallel framework where these conventional constraints are removed. Even in this case the systematics of the hadron yield ratios seems to completely overlap with the system including the constraints.
For the thermodynamic equilibration of the hadrons emerging from the heavy-ion collision experiments there could be many possibilities. For example, there could be the equilibration of different species of hadrons separately, or there could be the equilibration of groups of certain hadrons. The most interesting possibility is the equilibration of all the hadrons simultaneously under the umbrella of the conserved charges of strong interactions. The main intention of our present work was to find out the extent of such an equilibration achieved in the heavy-ion collision experiments. If the systematics studied here resulted in too large variations one would have to conclude that the holistic equilibration picture does not hold. However our results seem to indicate that for the full HRG spectrum chemical equilibration in terms of the conserved charges of strong interactions has been achieved to a very high degree for a wide range of √ s in the heavy-ion collisions.
Another outcome of our study is the variation of the NetQ/NetB ratio and net strangeness with √ s for the scheme without constraints. Surprisingly we find the strangeness neutrality to be still valid in the full range of √ s. This could only be achieved by an intricate adjustment of the strangeness chemical potential with respect to its counterparts for the other conserved charges. This reaffirms that all the three conserved charges work in unison to bring about the equilibration of hadronic matter in heavy-ion collisions.
We find the NetQ/NetB to remain at its intended value up to √ s ∼ 40 GeV, and thereafter generally decrease with √ s with growing uncertainties. This could be an indicator of the range of collision energies beyond which baryon stopping is reduced.

VIII. ACKNOWLEDGEMENTS
This work is funded by CSIR, UGC and DST of the Government of India. SB thanks A. Andronic and J. Stachel for their useful suggestions.