Chemical freeze-out systematics of thermal model analysis using hadron yield ratios

We provide a framework to estimate the systematic uncertainties in chemical freeze-out parameters extracted from $\chi^2$ analysis of thermal model, using hadron multiplicity ratios in relativistic heavy-ion collision experiments. Using a well known technique of graph theory, we construct all possible sets of independent ratios from available hadron yields and perform $\chi^2$ minimization on each set. We show that even for ten hadron yields, one obtains a large number ($10^8$) of independent sets which results in a distribution of extracted freeze-out parameters. We analyze these distributions and compare our results for chemical freeze-out parameters and associated systematic uncertainties with previous results available in the literature.


I. INTRODUCTION
Relativistic heavy-ion collisions provides the opportunity to create hot and dense QCD matter and study its thermodynamic and transport properties. The 'standard' model of relativistic heavy-ion collision has been developed in last few decades by analyzing the experimental data from Relativistic Heavy-Ion Collider (RHIC) at Brookhaven National Lab, USA and Large Hadron Collider (LHC) at CERN in Geneva, Switzerland. The analysis of experimental data suggests that the fireball produced in these collisions consists of deconfined quarks and gluons in the early stage of evolution which tends to thermalize rapidly. This quark-gluon plasma (QGP) undergoes a transition from partonic phase to hadronic phase during the later stage of evolution and finally the energy-momentum of these hadrons are measured by the detectors. The rapid thermalization of the fireball, due to strongly interacting constituents of the medium, provides the motivation to look for thermal description of the observed hadron yields. Thermal models of hadronic system has a long history [1][2][3][4][5]. Indeed, statistical models, with the assumption of complete thermalization of hadronic matter, has been quite successful in explaining the hadron yields and their ratios measured in relatively recent experiments .
Thermodynamic parameters obtained from statistical model analysis can be used to characterize the freeze-out hypersurface as the last surface of interaction. The onset of chemical freeze-out is said to have occurred when particle abundances are fixed and composition altering inelastic interactions have ceased. Subsequently kinetic freeze-out occurs when elastic interactions between particles have also stopped and the density drops to the level * sumana@jcbose.ac.in † a.jaiswal@niser.ac.in ‡ sutanu@niser.ac.in that the hadron momentum spectra remain unchanged.
Beyond kinetic freeze-out, hadrons are assumed to stream freely to the detector. Hadron resonance gas (HRG) model, which assumes a statistical description of mixture of hadrons and their resonances in thermodynamic equilibrium, leads to a good description of the medium at freeze-out, for a wide range of collision energy . In this model, thermodynamic equilibrium state of the strongly interacting matter is completely determined by the temperature (T ) and the three chemical potentials µ Q , µ B and µ S corresponding to baryon number (B), electric charge (Q) and strangeness (S), respectively. These parameters at freeze-out can be extracted from statistical model calculations by performing a χ 2 minimization fit to the available experimental multiplicity data [21][22][23][24][45][46][47][48][49][50][51]. Several codes like THERMUS [52], SHARE [53], THERMINATOR [54] are publicly available to compute the abundances of particles using such statistical hadronization approach.
The systematic uncertainties in the experimental data are expected to be reduced when one considers the ratios of hadron yields for χ 2 analysis, which overall cancels the the system volume (V ) [21,55]. However, there are multiple ways to form a set of (N − 1) independent ratios given N number of hadron yields. For simplicity, we call a set containing (N − 1) independent ratios as independent set. It is evident that the choice of particle ratios, may introduce a bias for the extracted freeze-out parameters [21,56]. Therefore, the freeze-out parameters extracted from a χ 2 minimization fit to these ratios would depend on the set under consideration. To avoid this problem, one may try to fit the absolute yields rather than their ratios which requires the inclusion of additional parameter V . This however is subject to its own bias [21,22,24]. It has also been observed that fitting absolute yield is more prone to converge to a false minima because of the strong correlation between yield normalization and other freeze-out parameters [53]. This problem becomes more pronounced when one tries to incorporate excluded volume effect [21]. On the other hand, yields ratios are not affected significantly when considering excluded volume effect and therefore the extracted freeze-out parameters remains stable.
In order to reliably extract freeze-out parameters from hadron yield ratios, one has to understand the systematic uncertainties due to the choice of specific ratios. First attempt towards understanding this uncertainty was made recently in Ref. [56] where the authors considered several independent sets for extraction of freeze-out parameters. These sets were chosen such that each of the N C 2 ratios appear in at least one set. However it is important to note that the choice of these particular sets are also not unique. Therefore, the bias arising from choice of those specific sets still remains which could only be removed by considering all possible independent sets. Given the importance of statistical models in the context of relativistic heavy-ion collisions, it is essential to conclusively address the issue of systematic uncertainty.
In this article, we provide a framework to estimate systematic uncertainties in the extraction of chemical freeze-out parameters from analysis of hadron multiplicity ratios. We construct all possible independent sets from available hadron yields by using a well known technique in graph theory. Subsequently, we perform χ 2 minimization on each set which leads to a distribution of the extracted freeze-out parameters. From these distributions, we obtain quantitative estimates of systematic uncertainty in the extracted freeze-out parameters corresponding to yield ratios of experimental data at 200 GeV (RHIC) and 2.76 TeV (LHC) collision energies. We also estimate these uncertainties after removing the usual constraints on the conserved charges. Finally, we compare our results for chemical freeze-out parameters and associated systematic uncertainties with previous results available in the literature.

II. HADRON RESONANCE GAS MODEL
We consider the HRG model for our analysis of yield ratios. The thermodynamic potential of HRG in terms of grand-canonical partition function is given by, where sum runs over all hadrons and resonances. The upper sign is for fermions and lower is for bosons, and the normalization factor V is the volume of the fireball. Here g i , E i and m i are respectively the degeneracy factor, energy and mass of i th hadron. While, µ i = B i µ B + S i µ S + Q i µ Q is the chemical potential of the i th hadron with B i , S i and Q i denoting its baryon number, strangeness and electric charge. For a thermalized system the number density n i ob-tained from partition function is given as, (2) The rapidity density for i'th detected hadron to the corresponding number density in the HRG model can be written as [23], where the subscript 'Det' denotes the detected hadrons.
If the heavier resonances j decay to the i-th hadron, then where n i and n j depend on the thermal parameters (T, µ B , µ Q , µ S ) which can be obtained by fitting experimental hadron yields to the model calculations. The fitted values of thermal parameters are obtained by minimizing the χ 2 , which is defined as, where, R exp α is the ratio of hadron yields obtained from experiments, R th α is the ratio of the number densities from the model calculations and σ exp α is the experimental uncertainty.
As individual multiplicities or their ratios are not conserved by strong interactions, χ 2 fitting is performed by choosing a set of multiplicity ratios. In general, with N given experimental hadron yields p 1 , p 2 , · · · , p N one can construct N (N − 1) ratios of the form p i /p j , where 1 ≤ i, j ≤ N and i = j. However, interchanging numerator and denominator does not provide any new information in our analysis. Therefore, it is sufficient to consider N (N −1)/2 ratios of the given particle yields. Because of the correlation existed between all possible multiplicity ratios, it is only relevant to choose (N − 1) number of statistically independent ratios out of those N (N − 1)/2 ratios to parameterize the thermal model [21,57]. It has also been pointed out within a thermal model analysis that, choosing a specific independent set may introduce bias in the minimization process [56]. However, it is possible to quantify this systematic uncertainty in the freeze-out parameters, extracted from χ 2 fitting, by considering all possible sets of independent ratios. Observe that, (N −1) ratios can be chosen in N (N −1)/2 C N −1 ways. But, any set containing (N − 1) ratios may not necessarily be independent. Also, enumerating all independent sets from N (N −1)/2 C N −1 sets, is a complex and tedious process because the total number of independent sets grows immensely for large N . We perform this enumeration process by relating it to a particular case of the minimum spanning tree enumeration problem, one of the well-known problems in combinatorial optimization.

III. GENERATING ALL INDEPENDENT SETS
Let S be a set of N distinct particle yields, namely p 1 , p 2 , · · · , p N . We associate an undirected graph G to S whose vertices are labeled by p 1 , p 2 , · · · , p N and there is an edge p i p j between the vertices p i , p j , if either p i /p j or p j /p i is in S. Observe that if p i p j , p j p k are edges of G then the vertices p i and p k are connected through the vertex p j . More generally, a path between two vertices p i and p j is said to exists if either there is an edge between them or they can be connected through some or all of the remaining N − 2 vertices. Since each of the particle yields p 1 , p 2 , · · · , p N must appear at least in one of the ratios in S, then every vertex of G is connected to at least one of the remaining N − 1 vertices. Let us denote the ratio of yields in S by r 1 , r 2 , · · · , r N −1 . Consequently, independence of S also implies that for each 1 ≤ k ≤ N − 1 The above equation implies that if p i p j is an edge of G then there is no other path between p i and p j . Equivalently, between any two vertices of G there is an unique path between them. Hence, G is a tree on the N vertices labeled by p 1 , p 2 , · · · , p N . Conversely, every tree on N vertices p 1 , p 2 , · · · , p N gives rise to an independent set. Therefore, we have a one to one correspondence between the trees on N vertices labeled by p 1 , p 2 , · · · , p N and the independent sets. For example, let us consider four distinct particles p 1 , p 2 , p 3 , p 4 . Then we have the following 16 independent sets: These sets can be obtained by associating a unique tree G j to every independent set I j for all 1 ≤ j ≤ 16 as shown in Fig. 1. Here the vertices symbolize particle yields. Each diagram in Fig. 1 corresponds to one set of independent ratios where the ratios are represented by lines joining two vertices. There are six distinct ratios (unique upto inverse) for N = 4 which are denoted by different colors. In general, the total number of different trees on N labeled vertices is equal to N N −2 , which is given by the Cayley's formula [58,59]. In fact the generation of those N N −2 many trees is equivalent to generating all spanning trees of the complete graph on N vertices labeled by p 1 , p 2 , · · · , p N . This is a particular case of the minimum spanning tree enumeration problem. This problem is important in its own right due to its wide range of applications including telecommunication networks, and therefore several algorithms with improvised efficiencies have been developed since mid 50's. We refer [60][61][62] and the references therein for more details.
The algorithm we have employed to generate all trees on N labeled vertices is as follows. It is straightforward to check that there can be N (N −1) 2 different edges connecting any two vertices p i , p j considering 1 ≤ i, j ≤ N and i = j. By implementing the algorithm, at every step we choose N − 1 edges say e 1 , e 2 , · · · , e N −1 (without repetition) from N (N −1) 2 edges and denote the graph formed with N vertices labeled by p 1 , p 2 , · · · , p N and edges e 1 , e 2 , · · · , e N −1 , by G ′ . In the next step, we shall verify whether the vertex p i , for all 1 ≤ i ≤ N , is connected to any other vertex by an edge or not. Then we construct the incidence matrix M = (m ij ) N ×(N −1) to verify whether the graph G ′ is a tree. The entries of M are defined as follows: m ij = 1 whenever the edge e j joins the vertex p i to some other vertex, and 0 otherwise, for all 1 ≤ i, j ≤ N . Observe that, if we delete one of the rows of M then we obtain a submatrix of M of size (N − 1) × (N − 1). We continue this process for each of the N rows of M and obtain N many submatrices of M . Then we check whether the rank of each of those N submatrices of M is N − 1. If it is true then we conclude that G ′ is a tree. In this way, we generate all N N −2 trees on N vertices and consequently, we obtain all N N −2 possible independent sets of ratios for our analysis.

IV. DATA ANALYSIS
We have used RHIC [63][64][65][66][67][68][69][70][71][72][73][74][75][76][77][78] and LHC [79][80][81][82] data at mid-rapidity for the most central collisions for our analysis. All hadrons with masses up to 2 GeV, with known degrees of freedom, are included for the HRG spectrum. The masses and branching ratios used are as given in Refs. [52,83]. It is to be noted that we have used vacuum masses for all the hadrons as the in-medium mass modification in connection with chiral symmetry restoration does not have any significant effect on freeze-out parameters for higher collision energies [84][85][86]. We could not find theΛ yield at LHC. Therefore at this energy, we assumedΛ yield to be same as that reported for Λ. Therefore, in order to obtain the freeze-out parameters, we have used the hadron yields of π ± (139. To obtain the freeze-out parameters it is a common approach to fix µ Q and µ S using these following constraint relations [45] , and With these two constraint equations, the problem is reduced to a two dimensional problem. The constant value of the ratio of net baryon number and net charge, depends on the physical system. For example, in Au + Au collisions, this constant is given by (N p + N n )/N p = 2.5, with N p and N n denoting the number of protons and neutrons in the colliding nuclei. However for LHC ( √ s N N = 2.76 TeV) and RHIC ( √ s N N = 200 GeV) energies, baryon stopping is expected to be negligible. In such cases, net baryon number and net charge may be vanishingly small and therefore their ratio may become arbitrary in a rapidity bin. Thus, it seems natural to relax the constraint given in Eq. (7) to extract freeze-out parameters χ 2 minimization. In the present analysis, we consider both cases: we call Model-I where the constraint given in Eq. (7) is imposed and Model-II where we do not enforce this constraint relation to hold. For N = 10 observed hadron yields, we have first obtained 10 8 independent sets of nine ratios, as explained in Sec. III. We then performed χ 2 minimization fits for each set using our numerical code for both Model-I and Model-II with a convergence criteria of 10 −4 or better. This procedure to fit all possible independent sets results in 10 8 values of each extracted freeze-out parameters and corresponding χ 2 numbers. We then construct histograms for each extracted freeze-out parameters, as well as for the χ 2 values. However, histograms give an overall notion about the density of the underlying distribution of the data. Therefore, for further analysis, we perform a least square fit to histograms of each parameter with Gaussian distribution. The fitted distribution curves quantify the probability distribution function for a population that has the maximum likelihood of producing the distribution that exists in the given sample. The mean of the Gaussian distribution corresponding to a given freeze-out parameter leads to an estimate of the central value of that parameter. The systematic uncertainties arising from the fitting procedure is calculated from full width at half maximum (2.3548σ), where σ is the standard deviation of the Gaussian distribution.

V. RESULTS AND DISCUSSIONS
In this section, we discuss the extracted chemical freeze-out parameters for center of mass collision energy √ s N N = 200 GeV at RHIC and √ s N N = 2.76 TeV at LHC. We fit 10 8 sets of yield ratios to obtain a distribution for fit parameters, i.e., freeze-out temperatures T and freeze-out chemical potentials µ B , µ Q and µ S . The estimated value of the parameters are derived from the mean of the fitted Gaussian distribution. The error bars quantifies the variations of the parameters obtained from all possible sets of independent ratios. Results including the constraint, Eq. (7), on net baryon and net charge (Model-I) and those for the unconstrained system (Model-II), are shown separately. Finally we also compare our results with those obtained earlier in literature.
Histogram distribution and corresponding Gaussian fit to this distribution for freeze-out temperatures are shown for center of mass energy  Fig. 2b. The qualitative as well as the quantitative estimations, for both Model-I and Model-II, are consistent with each other. We find that for RHIC the extracted mean value of freeze-out temperature is 169 MeV which is slightly larger than the earlier estimates. On the other hand, for LHC, we find the fitted mean value to be 155 MeV. Nevertheless the estimated values lies in near proximity with the values reported in literature, where the analysis was done with a particular set of independent ratios as documented in Table I and  in Table II. Our analysis shows that the present model has a small systematic uncertainty in chemical freeze-out temperature, which is quite promising. In general, the freeze-out temperature increases with increasing √ s N N and approaches a saturation [87], except at the LHC energy. As expected, in our analysis also, the temperature is lower at LHC energy than RHIC energy, due to the lower yield of protons at LHC [82].
In Fig. 3, we show the histogram distribution and the Gaussian fit to this distribution for baryon chemical potential µ B extracted from all independent sets. In Fig. 3a we present our results for √ s N N = 200 GeV collision energy whereas in Fig. 3b, results for √ s N N = 2.76 TeV is shown. Generally, for lower collision energies, significant number of baryons deposit their energies in the vicinity of center of mass frame due to baryon stopping. But as the energy increases colliding baryons may become almost transparent to each other and deposit their energy outside the collision region. In our analysis for both the models, we have estimated the value of µ B and quantifies the systematic uncertainty for high collision energies. dominant. The non-zero value of baryon chemical potential induces the production of strange baryons, which imminently requires the strange chemical potential to produce enough strange anti-particles. The existence of this non-zero strange chemical potential µ S , ensures the vanishing of net strangeness condition of the colliding nuclei. The appearance of non-zero µ S is considered as an indication of chemical equilibration in strongly interacting matter [88]. As the collision energy increases, µ S eventually decreases. On the other hand, electric charge chemical potential µ Q appears due to the net-isospin of the colliding nuclei. While scanning the collision energy range, µ Q remains non-zero but quite small and eventually approaches to zero at higher √ s N N . Several analyses thus approximated µ Q to be zero at higher collision energies [55]. In our analysis, we have estimated the value of µ S and µ Q using both Model-I and Model-II.
The mean value of µ S 's are in reasonable agreement with each other for both the models, as illustrated in Figs. 4a and 4b for RHIC energy. Though the estimated values commensurate for both the models, Model-II has much larger systematic uncertainty and exclusion of constraint is the source of this uncertainty. The same features are also observed for LHC energy as shown in Figs. 4c and 4d. Furthermore, as shown in Figs. 5a and 5b, the numerical value of µ Q is small. Especially for LHC energy, as shown in Figs. 5c and 5d, it is vanishingly small. Therefore, exclusion of net baryon to net charge constraint in Model-II result in fluctuations. Only the χ 2 analysis using these parameters can confirm the reliability of Model-II, which will be discussed later.
For comparison, the freeze-out parameters from relevant literature are also documented in Tables I and II  for collision    must keep in mind that there are significant differences in the analysis procedure and computation techniques employed here. The error bars considered in Ref. [21] are obtained by varying the χ 2 by unity, which accounts for the fitting error. In Ref. [55], µ Q was assumed to be vanishingly small and therefore neglected. However, an additional parameter γ s was introduced to account for the possible deviation of strange particle abundances from chemical equilibrium. The only set of ratios that are used in the analysis of Ref. [55] are π − /π + , k − /k+, p/p,Λ/Λ, Ξ + /Ξ − , k − /π − ,p/π − , Λ/π − and Ξ + /π − . They have also studied different sources of systematic uncertainties, which propagates from yield to freeze-out parameters. In Ref. [56] the systematic error bars has been calculated with nine independent sets, considering the usual constraint of strong charges (WC) and without constrain relations (WOC) of Eq. (7) and Eq. (8).
In Ref. [40], the authors have considered system volume as a parameter and the standard deviations quoted there comes only from experimental uncertainties. Reduced χ 2 i.e. χ 2 /NDF gives a measure of goodness of fit, where NDF stands for the number of degrees of freedom. The usual practice is to choose a specific independent set so that the χ 2 is equal to the number of degrees of freedom of the system. The histogram distribution of reduced χ 2 corresponding to all possible independent sets for √ s N N = 2.76 TeV (LHC) and √ s N N = 200 GeV (RHIC) are shown for both Model-I and Model-II in Fig. 6. We see that χ 2 /NDF for both the models differ only slightly. It is interesting to see that the distribution is peaked at χ 2 /NDF ≃ 1 for RHIC as shown in Fig. 6a indicating a good fit. On the other hand, from Fig. 6b, we see that χ 2 /NDF ≃ 2.2 for LHC. This suggests that quality of fit at LHC may not be as good as that at RHIC. Anyhow, this result is not surprising as it has been pointed out in existing literature that reduced χ 2 for LHC is quite large, if we consider the equilibration of all the hadrons simultaneously within the framework of the conserved charges [46,89]. However, the usual practice is to minimize χ 2 for a specific independent set and large value of χ 2 /NDF for that particular set does not necessarily be alarming. But, we have done the study for all possible independent sets and observe that the peak of the reduced χ 2 distribution lies near 2.2, for both the models. Thus, the thermal model analysis for LHC data demands further study.

VI. SUMMARY AND OUTLOOK
It has been conjectured in existing literature that χ 2 fit of experimentally measured hadronic yields, or multiplicity ratios, within thermal model framework can successfully extract chemical freeze-out parameters to a certain extent. It is evident that, there are many possible sources of systematic uncertainties that can affect the analysis. The usual practice is to choose a specific set of independent ratios for the analysis, such that χ 2 /NDF is close to one. This is one of the sources of systematic uncertainty as the equilibration parameters are biased to the chosen particle ratios. Estimation of these uncertainties are important because it would shed light on the reliability of such thermal models where we have assumed a grand canonical ensemble with thermodynamic equilibration of all the hadrons emerging from the heavy-ion collision experiments.
In this article, we have provided an elegant method to quantify such systematic uncertainties in the chemical freeze-out parameters extracted from thermal model analysis of hadron multiplicity ratios. We have identified the enumeration problem of independent sets to a well known problem in combinatorial optimization and graph theory. By implementing the minimum spanning tree algorithm, we have generated all possible independent sets. This corresponds to the number of all different trees on N vertices resulting in generation of N N −2 independent sets as expected from Cayleys formula.
We have performed the χ 2 minimization on each sets to obtain freeze-out parameters corresponding to yield ratios of experimental data at collision energies of √ s N N = 200 GeV at RHIC and √ s N N = 2.76 TeV at LHC. Since the number of sets are extremely large (10 8 ), we obtain a distribution of these parameters. From these distribu-tions, we estimate the mean value of the said parameter along with the quantitative measure of systematic uncertainty. Model-II seems to have larger systematic uncertainty than Model-I, which is expected because Model-II has one extra free parameter. But the χ 2 /NDF distribution for RHIC is peaked around one which ensures that the systematics are under control for both the models. One may be tempted to conclude that a global equilibration scenario is achieved for RHIC where equilibration of all hadrons occurs simultaneously. On the other hand, for LHC, we found that although the parameters seem to be in good agreement with existing literature, the value where of χ 2 /NDF distribution is peaked is away from one. This leads to a serious challenge about the applicability of the thermal model (in the present form) at LHC energy.
Looking forward, the framework presented here can be applied to other collision energies at RHIC and LHC in order to ascertain that the extraction of chemical freezeout parameters are not biased by the choice of yield ratios. Moreover, the current framework can be easily extended to physics analysis using double/multiple ratios which will have important implications for sequential freeze-out models. We leave these for future work.