On the Higher Moments of Particle Multiplicity, Chemical Freeze-Out and QCD Critical Endpoint

We calculate the first six non-normalized moments of particle multiplicity within the framework of the hadron resonance gas model. In terms of the lower order moments and corresponding correlation functions, general expressions of higher order moments are derived. Thermal evolution of the first four normalized moments and their products (ratios) are studied at different chemical potentials {\mu}, so that it is possible to evaluate them at chemical freeze out curve. It is found that a non-monotonic behavior, reflecting the dynamical fluctuation and strong correlation of particles starts to appear from the normalized third order moment. We introduce novel conditions for describing the chemical freeze out curve. Although the hadron resonance gas model does not contain any information on the criticality related to the chiral dynamics and singularity in some physical observables, we are able find out the location of the QCD critical endpoint at $\mu\sim 350 $MeV and temperature $T \sim 162 $MeV.


I. INTRODUCTION
Recently, the higher order multiplicity moments have gained prominence in high energy physics with a huge hope in pinpointing the QCD critical endpoint (CEP) [1] connecting the first order boundary separating the hadronic from the partonic matter at high density with the crossover boundary at low density [2,3]. It is clear that the hadron resonance gas partition function is an approximation to a non-singular part of the free energy of QCD in the hadronic phase. There are large theoretical uncertainties of its location and even its existence is not fully confirmed, yet [4][5][6]. The characteristic feature of CEP are critical dynamical fluctuations [7][8][9][10][11]. The higher order moments are conjectured to reflect the large fluctuations associating the hadronquark phase transition. This was the motivation of a remarkable number of experimental and theoretical studies [1, [12][13][14]. Recently, various calculations have shown that the higher order moments of the multiplicity distributions of some conserved quantities, such as net-baryon, net-charge, and net-strangeness, are sensitive to the correlation length ξ [15][16][17], which in turn is related to the higher order moments, themselves. In realistic heavy-ion collisions, the correlation length is found to remain finite.
Apparently, it should not be an exception that the strongly interacting QCD matter undergoes phase transition(s) at extreme conditions, as almost all matter types suffer from such a critical change as the extreme conditions change. The Lattice QCD calculations predict that a cross-over takes place between the hadronic phase and the Quark Gluon Plasma (QGP) phase, when the temperature exceeds critical value of T c ≃ 150 − 190 MeV. As given in Ref. [3,18], depending on the different parameters (for instance, the quark favours and their masses) and on the lattice configurations, lattice QCD assigns different values to T c . Furthermore, the value of T c depends on the baryon chemical potential µ. With vanishing µ the lattice QCD calculations [19] show that the higher order susceptibilities can be related to the higher order moments of the corresponding multiplicity distributions. They also show non-monotonic behavior near T c .
Apparently, the QCD phase diagram including CEP can be mapped out using this characteristic behaviour. In literature, one finds that the QCD phase diagram has been determined using various criteria, such as critical energy density [20][21][22] or line of constant physics [3]. Also the chemical freeze-out curve can be characterized assuming constant physics. For a recent and extensive review, we refer to [23] and the references therein.
At large chemical potential µ, various calculations based on QCD-based models indicate that the transition from the hadronic phase to the QGP phase is of first order. The endpoint connecting this line with the one of the cross over (at small µ) is likely of second order [4-6, 6, 25]. Therefore, the phase transition should be accompanied with large fluctuations in different physical quantities. Experimentally, we study the QCD phase diagram by varying the colliding energy in heavy-ion collisions and determining the critical temperature T c [2]. The possibility of finding QCD CEP has motivated the RHIC beam energy scan program [31]. By tuning the collision energy from a center-of-mass energy of 200 GeV down to 5 GeV, we are able to vary the baryon chemical potential from ∼ 2 to ∼ 500 MeV. Nevertheless, there are large theoretical uncertainties of its location and even its existence is not fully confirmed [4][5][6]24]. The higher order moments of various physical quantities are conjectured to reflect the non-monotonic behaviour near the QCD CEP. In the present work, we give a systematic study for the higher order moments of the multiplicity distribution and show their ability to represent the fluctuations along the QCD phase diagram.
Many reasons speak for utilizing the physical resonance gas model (HRG) in predicting the hadron abundances and their thermodynamics. The HRG model seems to provide a good description for the thermal evolution of the thermodynamic quantities in the hadronic matter [3, 7-9, 20-22, 26, 27] and has been successfully utilized to characterize the conditions deriving the chemical freeze-out at finite densities [28,29]. In light of this, HRG can be used in calculating the higher order moments of particle multiplicity using a grand canonical partition function of an ideal gas with all experimentally observed states up to a certain large mass as constituents. The HRG grand canonical ensemble includes two important features [3]; the kinetic energies and the summation over all degrees of freedom and energies of resonances.
On other hand, it is known that the formation of resonances can only be achieved through strong interactions [30]; Resonances (fireballs) are composed of further resonances (fireballs), which in turn consist of resonances (fireballs) and so on. In other words, the contributions of the hadron resonances to the partition function are the same as that of free particles with some effective mass. At temperatures comparable to the resonance half-width, the effective mass approaches the physical one [3]. Thus, at high temperatures, the strong interactions are conjectured to be taken into consideration through including heavy resonances. It is found that hadron resonances with masses up to 2 GeV are representing suitable constituents for the partition function [3, 7-9, 20-22, 26, 27]. Such a way, the singularity expected at the Hagedorn temperature [20,21] can be avoided and the strong interactions are assumed to be considered.
Nevertheless, the validity of HRG is limited to temperatures below the critical one, T c .
In this paper, we study the first six non-normalized higher order moments of the particle multiplicity in the hadron resonance gas (HRG) model, which is introduce in section II. In section II A, we suggest an answer to the question: What is needed when going from "trivial" lower order to "sophisticated" higher order moments? Based on results, general expressions for arbitrary higher moments are deduced, so that it is possible to conclude that going from lower to higher order moments is achievable through a series of lower order moments and correlation functions. The thermal evolution of the first four normalized moments and their products (ratios) are studied at different chemical potentials µ in section II B. Therefore, it is possible to evaluate them at the chemical freeze-out curve, which is characterized by constant s/T 3 at all values of µ, where s and T are entropy density and temperature, respectively. Sections III and IV are devoted to introduce novel conditions describing the chemical freeze-out curve and define the location of the QCD critical endpoint.

II. HIGHER ORDER MOMENTS OF THE PARTICLE MULTIPLICITY
In a grand canonical ensemble, it is straightforward to derive an expression for the pressure.
In section II A, we give a list of different moments of particle multiplicity. As we move from lower order to higher order multiplicities, certain distribution functions are being added / subtracted. These distribution functions represent higher order correlations. In section II B, we list out various normalized higher order moments, the normalization being based on the variance σ. The products (ratios) of normalized moments are elaborated in section II C.
The hadron resonances treated as a free gas [3,[20][21][22]26] are conjectured to add to the thermodynamic pressure in the hadronic phase (below T c ). This statement is valid for free as well as for strongly interacting resonances. It has been shown that the thermodynamics of strongly interacting system can also be approximated to an ideal gas composed of hadron resonances with masses ≤ 2 GeV [3,32]. Therefore, the confined phase of QCD, the hadronic phase, is modelled as a non-interacting gas of resonances. The grand canonical partition function reads where H is the Hamiltonian of the system and T is the temperature. The Hamiltonian is given by the sum of the kinetic energies of relativistic Fermi and Bose particles. The main motivation of using this Hamiltonian is that it contains all relevant degrees of freedom of confined and strongly interacting matter. It includes implicitly the interactions that result in resonance formation. In addition, it has been shown that this model can submit a quite satisfactory description of particle production in heavy-ion collisions. With the above assumptions the dynamics the partition function can be calculated exactly and be expressed as a sum over single-particle partition functions Z 1 i of all hadrons and their resonances.
where ǫ i (k) = (k 2 + m 2 i ) 1/2 is the i−th particle dispersion relation, g i is spin-isospin degeneracy factor and ± stands for bosons and fermions, respectively.
Before the discovery of QCD, a probable phase transition of a massless pion gas to a new phase of matter was speculated [33]. Based on statistical models like Hagedorn [34] and Bootstrap [35], the thermodynamics of such an ideal pion gas is studied, extensively. After the QCD, the new phase of matter is known as quark gluon plasma (QGP). The physical picture was that at T c the additional degrees of freedom carried by QGP are to be released resulting in an increase in the thermodynamic quantities like energy and pressure densities. The success of hadron resonance gas model in reproducing lattice QCD results at various quark flavours and masses (below T c ) changed this physical picture drastically. Instead of releasing additional degrees of freedom at T > T c , the interacting system reduces its effective degrees of freedom at T < T c . In other word, the hadron gas has much more degrees of freedom than QGP.
At finite temperature T and baryon chemical potential µ i , the pressure of the i-th hadron or resonance species reads As no phase transition is conjectured in HRG, summing over all hadron resonances results in the final thermodynamic pressure in the hadronic phase.
The switching between hadron and quark chemistry is given by the relations between the hadronic chemical potentials and the quark constituents; µ i = 3 n b µ q + n s µ S , where n b (n s ) being baryon (strange) quantum number. The chemical potential assigned to the light quarks is µ q = (µ u + µ d )/2 and the one assigned to strange quark reads µ S = µ q − µ s . The strangeness chemical potential µ S is calculated as a function of T and µ i under the assumption that the overall strange quantum number has to remain conserved in heavy-ion collisions [3].
A. Non-normalized higher order moments As given in the introduction, the higher order moments can be studied in different physical quantities. For example, the higher order moments of charged-particle multiplicity distribution have been predicted four decades ago [36]. Recently, the higher order moments of various multiplicity distributions have been reported by STAR measurements [12,13] and lattice QCD calculations [37][38][39]. The empirical relevance of the higher order moments to the experimental measurements has been suggested in Ref.
[1]. Accordingly, the measurement of the correlation length seems to be very much crucial. In a future work, we shall discuss the experimental signatures of the higher order moments. Another ingredient would be the experimental sensitivity for the suggested signatures, which are based on singular behavior of thermodynamic functions. For the time being, we just want to mention that the experimental measurements apparently take place at the final state of the collision, which would mean that the signals have to survive the extreme conditions in such high energy collisions. It has been pointed out that the contribution of the critical fluctuations to the higher order moments is proportional to a positive power of ξ. The latter is conjectured to diverge at QCD CEP and such an assumption is valid in the thermodynamic limit.
In the present work, we first study the contribution of the hadron resonance mass spectrum in clarifying the interrelations between the various non-normalized and normalized moments.
The motivation of using resonance gas is clear [20][21][22]. As introduced about five decades ago, the resonance gas is only valid below T c [3, 20-22, 26, 34]. Therefore, the present work is valid in the hadronic phase, only.
For the i-th particle, the "first" order moment is given by the derivative of p = Thus, from Eqs. (4)-(9), a general expression for the r-th order moment can be deduced where the coefficients read where l ≤ r and a r,l vanishes ∀ r < 1. It is obvious that the coefficients of a certain moment are to be determined from a long chain of all previous ones. Such a conclusion dates back to about four decades [36], where it has been shown that the coefficients are related to high order correlation functions. Should this assumption is proven to be valid, then expression (11) gets a novel interpretation. It seems to sum up the correlation functions up to the r-th order.
According to [36] and when neglecting the two-particle correlations C 2 , then the higher order moments read where δN = N − N , and p i is the i−th particle. To the second order moment m 2 we have to add the effects of the two particle correlation function, 2 C 2 . The third order moment m 3 gets approximately three times this amount. The three and four particle correlations, Eqs. (15) and (16), appear first in third and fourth moment, Eq. (14) and (15), respectively.
When ignoring these high order correlations, we are left with the widely used two particle correlations. In the present work, we restrict the discussion to this type, only. It is of great interest, as it is accessible experimentally and achievable, numerically. The two particle correlations are suggested as a probe for the bulk QCD medium, energy loss, medium response, jet properties and intensity interferometry [41][42][43][44][45][46][47]. In addition to this list of literature, the comprehensive review [48] can be recommended. Taking into account the particle multiplicities, then expressions (12), (13) and (14) can be re-written as follows.
As mentioned in section II A, the experimental signals for the higher order moments shall be discussed in a future work.
As given in Eq. (10), the dependence of first six non-normalized moments on lower ones can explicitly be deduced as follows.
Naively spoken, we conclude that raising lower order moments to higher ones is achievable through a series of all lower order moments and an additional term reflecting the correlations, themselves [36]. From Eqs. (20)- (25), the additional terms are proportional to N r . Apparently, they are essential in order to judge whether going from lower to higher order moments would make the signatures of dynamical fluctuations clearer than when excluding them.
The first and second terms can be generalized. Then, a general expression would read where c mr is the coefficient of the r-th moment. The last term in Eq. (26) can be elaborated when higher order moments are calculated. The latter are essential in order to find out a clear pattern.
The figure 1 presents the results of the higher order moments as calculated in the HRG model.
As introduced previously, raising the orders of multiplicity moments results in new coefficients and new integrals. The earlier are partly characterized in Eq. (26). The dimensionless nonnormalized moments are given in dependence on the temperature T in the left panel of Fig.   1. The new coefficients that resulted in when going from lower to higher order moments are depicted in the right panel. The integrals can be compared with the correlation functions [36].
We find that these integrals are related to N r , where r is the order of the moment. It is obvious that successive moments have a difference of about one order of magnitude. Therefore, higher order N r can be disregarded with reference to their vanishing contributions.
In order to draw a picture about the contributions of the integrals appearing in Eqs. (20)- (25), the phase space integral can be replaced by a series representation. For the i-th particle species, the pressure can be expressed as where K 2 is modified Bessel function of the second kind and the factor (±) n+1 represents fermions and bosons, respectively. In deriving previous expression, we expand the functions appearing in Eq. (4), for instance. The resulting expression is applicable as long as µ i < m i is valid. Then, the derivatives with respect to µ/T leads to which is valid for all order moments of order r ≥ 1.
For completeness, we mention here that taking into consideration one particle and its antiparticle is not rare in literature. In this case, the expression (27) can be re-written as and accordingly, Comparing Eq. (28) with Eq. (27), we conclude that calculating higher order moments is associated with a gradual change in n. Expression (10) reflects almost the same behaviour, where higher order moments are associated with change in the coefficients.

B. Normalized higher order moments
The normalization of higher order moments which can be deduced through derivatives with respect to the chemical potential µ of given charges, apparently gives additional insights about the properties of higher order moments. From statistical point of view, the normalization is done with respect to the standard deviation σ, which is be related to ξ. Therefore, it provides with a tool to relate moments with various orders to the experimental measurement. The susceptibility of the distribution give a measure for σ. Should we are interested on the multiplicity, then the susceptibility is simply given by the derivative of first order moments with respect to µ. It has been shown that the susceptibility is related to ∼ ξ 2 [1]. The results of σ in hadronic, bosonic and fermionic resonances are calculated at different µ and given in Fig. 2. As discussed earlier, the strange quantum number has to remain conserved in high energy collisions. As per the standard model, this is one of the global symmetries in strong interactions. The procedure of keeping strange degrees of freedon conserved in HRG is introduced in Ref. [3]. This is the origin of the µ-dependence, especially in the bosonic gas. Although, baryon chemical potential µ vanishes per definition, the chemical potential associated with strange quark µ S remains finite. The dependence of µ S on µ is depicted in Fig. 16. Another feature in these calculations is the assumption that the freeze-out boundary is determined by constant s/T 3 , where s is the entropy density [28]. At chemical freeze-out boundary, the dependence of σ on µ is given in  The "third" order moment normalized to σ 3 is known as 'skewness". For standard Gaussian distribution, the skewness is obviously vanishing. Therefore, the skewness is an ideal quantity probing the non-Gaussian fluctuation feature as expected near T c . The QCD CEP is conjectured to be sensitive to skewness. Experimentally, it has been shown that the skewness S is related to ∼ ξ 4.5 [15]. The skewness for bosonic and fermionic resonance gas, respectively, reads At different chemical potentials, the skewness S is calculated in dependence on T and given in  The normalization of 4−th order moment is known as heteroskedacity or kurtosis. It means varying volatility or more accurately, varying variance. Actually, the kurtosis is given by normalized 4−th order moment minus 3. The subtraction of 3, which arises from the Gaussian distribution, is usually omitted [40,49,50]. Therefore, the kurtosis is an ideal quantity for probing the non-Gaussian fluctuation feature as expected near T c and CEP. A sign change of skewness or kurtosis is conjectured to indicate that the system crosses the phase boundary [19,51]. As HRG is valid below T c , the sign change is not accessible. It has been shown that kurtosis κ is related to ∼ ξ 7 [15].
C. Products of higher order moments There are several techniques to scale the correlation functions. The survey system's optional statistics module represents the most common technique i.e., Pearson or product moment correlation. This module includes the so-called partial correlation which seems to be useful when the relationship between two variables is to be highlighted, while effect of one or two other variables can be removed. In the present work, we study the products of higher order moments of the distributions of conserved quantities. The justification of this step is that certain products can be directly connected to the corresponding susceptibilities in Lattice QCD simulation and related to long range correlations [19,52,53]. Seeking for simplicity, we start with the Boltzmann approximation.

Modified Bessel function: Boltzmann statistics
When relativistic momentum integrals are replaced by summation over modified Bessel functions, the products of higher order moments in Boltzmann approximation read The origin of the second term in all expressions containing κ is obvious. Expressions (35)- (38) justify the conclusions in [54] that in Boltzmann approximation while These expressions are valid in the final state, which can be characterized by chemical and thermal freeze-out. In other words, these expressions are functions of the chemical potential.
Using relativistic momentum integrals (see section II C 2) shows that the products of moments result in a constant dependence on T [54].

Relativistic Momentum Integrals; quantum statistics
The fluctuations of conserved quantities are assumed to be sensitive to the structure of the hadronic system in its final state. As mentioned above, crossing the phase boundary or passing through critical endpoint is associated with large fluctuations. Most proposed fluctuation of observables are variations of second order moments of the distribution, such as particle ratio [11] and charged dynamical measurement [55]. Then, the fluctuations are approximately related to ξ 2 [17].  The ratio of standard deviation σ 2 and the mean multiplicity N for fermions and bosons reads where ± stands for fermions and bosons, respectively. The results are given in Fig. 4. We notice that the bosonic resonance gas results in smaller values than the fermionic one, especially in final state. Furthermore, we notice that σ 2 / N decreases with increasing µ of the bosons at middle temperature. This is exactly the opposite in the fermionic resonance gas. This would give an explanation for the observation that the results in hadron resonances are not spread as in the other two sectors. Figure 5 shows the dependence of σ 2 / N on µ at the chemical freeze-out boundary, which is characterized by s/T 3 = 7. The ratio σ 2 / N is equivalent to m 2 /m 1 , Eq. (48). The dependence of bosonic and fermionic σ on µ is given in Fig. 5, as well.
We notice that σ is smaller than σ 2 / N , especially at small µ. At large µ both quantities are almost equal.   The multiplication of skewness S by the standard deviation σ is directly related to the thermodynamics of the number susceptibility of the lattice QCD. In HRG, the bosonic and fermionic products read The product S σ is equivalent to m 3 /m 2 , Eq. (49). The results are given in Fig. 6. It is obvious that S σ ≃ 1 for either bosons or fermions. Then, for hadrons, S σ ≃ 2. Nevertheless, the fine structure seems to reveal interesting features.
The dependence of σ and σ 2 / N on T is illustrated in Figs. 2 and 4, respectively. It is obvious that both quantities have a monotonic behavior. Their dependences on µ are given in Fig. 5. Also, this type of dependences seems to be monotonic. As given in Fig. 6, the product S σ has a characteristic dependence on T . Regardless, the tiny change, we notice that increasing µ increases the bosonic S σ product, but decreases the fermionic S σ product. In both cases, it forms a folding fan. The hadronic product makes an amazing bundle in the middle (at a characteristic T ). This will be discussed in details in section III. When S is studied as a function of µ and given in Fig. 15, we find that S(µ) b remains almost constant, while S(µ) f raises with increasing µ. We also notice that both curves cross at a certain point. This is not the case of the µ-dependence of σ and σ 2 / N , Fig. 5. So far, we conclude that starting from the third normalized moment, a non-monotonic behavior appears. Thus, the kurtosis and its products should be expected to highlight such a non-monotonic behavior [19,52,56].
The multiplication of kurtosis by σ 2 called κ ef f [57] is apparently equivalent to the ratio of 3-rd order moment to 2-nd order moment, Eq. (50). In lattice QCD and QCD-like models, κ ef f is found to diverge near the critical endpoint [19,52]. In HRG, the bosonic and fermionic products read When ignoring the constant term in Eqs. (33) and (34), then the second terms in the previous expressions entirely disappear. The results are given in Fig. 7. In the hadronic sector, the dependence of κ σ 2 on the temperature T at different µ-values is given in the left panel. We notice that increasing T is accompanied with a drastic declination in κ σ 2 . Also, we find that κ σ 2 flips its sign at large T . When comparing the thermal evolution of κ σ 2 with the one of S σ, Fig. 6, we simply find that the latter is much drastically changing than the earlier. Also, when comparing their dependences on the chemical potentials at the freeze-out boundary, Figs. 9 and 14, it is apparent that the µ-dependence increases when the normalized fourth order moment is included. The product κ σ 2 calculated at the freeze-out boundary leads to some interesting findings. First, κ σ 2 almost vanishes or even flips its sign. Second, the T and µ corresponding to vanishing κ σ 2 are coincident with the phenomenologically measured freeze-out parameters.
Third, the freeze-out boundaries of bosons and fermions are crossing at a point located very near to the one assumed by the lattice QCD calculations to be the QCD CEP, section IV. The skewness and kurtosis can be combined through the standard deviation as follows.
Also when ignoring the constant term in Eqs. (33) and (34), the second terms in previous expressions disappear. The results of κ σ/S are given in Fig. 8. In the fermionic sector, κ σ/S decreases with increasing T . Increasing µ makes the decrease much faster. Although the drastic change in κ σ/S, this behavior can be compared with S σ, Fig. 6, qualitatively. In the bosonic sector, there is a very slow decrease with increasing µ. The same sector in S σ, Fig. 6, shows a very slow but an increasing dependence on T . Therefore, in the hadronic sector, κ σ/S is overall diminishing with increasing T . The declination becomes faster for larger µ. Fig. 9 illustrates the µ-dependence, i.e. κ σ/S is estimated at the freeze-out curve. We notice that the bosons result in an almost unvarying κ σ/S with growing µ. Therefore, the overall decrease in the hadronic sector is originated in the fermionic degrees of freedom.
Obviously, other products would complete missing ratios, for example, relating third and fourth order moments to the first one, respectively. As discussed above, the susceptibility χ is equivalent to σ 2 . These products seem to be sensitive to the volume independent multiplicities. The results of S χ ≡ S σ 2 are given in Fig. 10. At the freeze-out boundary, the µ-dependence of S σ 2 seems to be much stronger than the µ-dependence on S σ, Fig. 11.
The thermal evolution of κ σ 4 / N is illustrated in Fig. 12. Increasing T is accompanied by an almost overall increase of κ σ 4 / N values. At low T , fluctuations appear. In bosonic sector, At very low µ, it turns out to small positive values.

III. CHEMICAL FREEZE-OUT
In rest frame of produced particle, the hadronic matter can be determined by constant degrees of freedom, for instance, s/T 3 (4/π 2 ) = const [58]. The chemical freeze-out is related to the particle creation. Therefore, the abundances of different particle species are controlled by chemical potential, which obviously depends on T . With the beam energy, T is increasing, while the baryon densities at mid-rapidity is decreasing. The estimation of the macroscopic parameters of the chemical freeze-out can be extracted from particle ratio. These parameters collected over the last two decades seem to fellow regular patterns as the beam energy increases [23,28]. The higher order moments have been suggested to control the chemical freeze-out, so that several conditions have been proposed [59].
As introduced in section II C 2, the thermal evolution of κ σ 2 declines. This result seems to update previous ones [38,60], where κ σ 2 is assumed to remain finite and positive with increasing µ. In the present work, we find that the sign of κ σ 2 is flipped at high T [59]. It has been found that the T -and µ-parameters, at which the sign is flipped are amazingly coincide with the ones of the chemical freeze-out, Fig. 14. Also, we find that the freeze-out boundaries of bosons and fermions are crossing at one point located at the hadronic curve. This point is close to the one that lattice QCD calculations suggest for QCD CEP, section IV. The results are given in Fig. 14. Vanishing κ σ 2 leads to The rhs and lhs in both expressions can be re-written as which is valid for bosons and fermions. Then, the chemical freeze-out is defined, if the condition is fulfilled. At chemical freeze-out curve, a naive estimation leads to ξ ∼ 3 1/3 fm. In doing this, it is assumed that the proportionality coefficients of κ ∼ ξ 7 and χ ∼ ξ 2 , are equal. In the heavy-ion collisions, ξ has been measured [61]. Near a critical point, the experimental value ∼ 2 − 3 fm (only factor 3 larger) agrees well with our estimation.
At the chemical freeze-out curve, the intensive parameters T and µ which are related to the extensive properties entropy and particle number, respectively, have to be determined over a wide range of beam energies. Left panel of Fig. 14 collects a large experimental data set. For a recent review, we refer to [23] and the references therein. The double-dotted curve represents a set of T and µ, at which κ σ 2 vanishes. It is obvious that this curve reproduces very well the experimental data. As given above, at this curve the normalized fourth order moment κ is equal to three times the squared susceptibility χ. This new condition seems to guarantee the condition introduced in [28]; s/T 3 = const., right panel of Fig. 14, at least over the range 0 < µ < 0.6 GeV. When excluding all degrees of freedom but the fermionic ones, s/T 3 decreases with increasing µ (dashed curve). An opposite dependence is accompanied with the bosonic degrees of freedom (dotted curve).

IV. QCD CRITICAL ENDPOINT
In right panel of Fig. 14, it is interesting to notice that both fermionic and bosonic curves intersect at the hadronic curve at one point. Also, in the freeze-out diagram, all curves crossing at one point, left panel of Fig. 14. Furthermore, based on non-perturbative convergence radius, the critical endpoint as calculated in lattice QCD [5,6] given by solid square in Fig. 14 is located very near to the crossing point. It is clear that the HRG model does not contain any information on criticality related with the chiral dynamics and singularity in physical observables required to locate the CP. The HRG partition function is an approximation to a non-singular part of the free energy of QCD in the hadronic phase. It can be used only as the reference for LGT calculations or HIC to verify the critical behavior, but it can not be used as an origin to search for the chiral critical structure in QCD medium. This is the motivation of Fig. 15. and/or dependently, in order to determine µ of the crossing point: Due to the mathematical difficulties in dealing with these expressions, the integral over phase space has to be simplified. A suitable simplification is given in section II C 1.
where n b and n s being baryon (strange) quantum number and µ q (µ s ) is the baryon (strange) chemical potential of light and strange quarks, respectively. The quarks chemistry is introduced in section II. Accordingly, the difference between baryons and fermions is originated in the exponential function. For simplicity, we consider one fermion and one boson particle. Then, the baryon chemical potential µ at the chemical freeze-out curve, at which the fermionic and bosonic skewness (or kurtosis) curves of these two particles cross with each other can be given as In the relativistic limit, K 2 (m/T ) ≈ 2T 2 /m 2 −1/2 while in the non-relativistic limit K 2 (m/T ) ≈ πT /2m exp(−m/T )(1 + 15T /8m). It is obvious that the bosonic and fermionic degrees of freedom play an essential role in determining Eq. (61). Furthermore, it seems that the chemical potential of strange quark has no effect at the crossing point.
The dependence of µ S on µ b as calculated in HRG is given in Fig. 16. As mentioned above, µ S is calculated to guarantee strange number conservation in heavy-ion collisions. At small µ b , µ S has a linear dependence, µ S = 0.25 µ b (Hooke's limit). At large µ b , the dependence is no longer linear. Hooke's Limit Fig. 16: µ S as a function of µ b in hadron resonance gas model (dotted curve). The linear fitting is given by solid curve.
As discussed earlier, the critical behavior and the existence of QCD CEP can be identified by means of signatures sensitive to singular parts of the free energy, especially the ones reflecting dynamical fluctuations of conserved charges, such as baryon number and charge density [62].
The reason that the fermionic and bosonic higher order moments are crossing on a point very near the QCD CEP shall be elaborated in a forthcoming work [63].

V. CONCLUSION
In the present work, the first six non-normalized order moments of the particle multiplicity are calculated in the hadron resonance gas model. A general expressions for arbitrary higher order moments are deduced in terms of lower ones. We concluded that going from lower to higher order moments is possible through adding up a series consists of alls lower order moments plus correlation functions. We studied the thermal evolution of the first four normalized order moments and their products (ratios) at different chemical potentials µ. By doing that, it was possible to evaluate first four normalized moments at the chemical freeze-out curve. The freeze-out curve is characterized by constant s/T 3 at all values of µ, where s and T are entropy density and temperature, respectively. It has been found that non-monotonic behavior reflecting dynamical fluctuation and strong correlations appears starting from the normalized third order moment (skewness S). Furthermore, non-monotonicity is observed in the normalized fourth order moment, the kurtosis κ, and its products. These are novel observations. Although HRG is exclusively applicable below T c , i.e. it does not include deconfinement phase transition, it is apparent that the higher order moments are able to give signatures for the critical phase transition. The accuracy of our calculations made it possible to have clear evidences on the phase transition. Based on these findings, we introduced novel conditions characterizing the chemical freeze-out curve and the location of the QCD critical endpoint as follows. The chemical freeze-out curve is described by κ = 3 χ 2 , where χ is the susceptibility in particle number, i.e. the second order moment. The location of QCD critical endpoint (at T -and µ-axis) is positioned when the condition S b = S f or κ b = κ f is fulfilled. The subscripts b and f refer to bosons and fermions, respectively. Accordingly, the QCD endpoint is positioned at µ ∼ 350 MeV and T ∼ 162 MeV. We are able to estimate these two quantities, although the hadron resonance gas model basically does not contain information on criticality related with the chiral dynamics and singularity in physical observables required to locate the critical endpoint. After submitting this work, a new paper posted in arXiv, in which the authors used the second order moment (susceptibility) and the fourth order one (kurtosis) are apparently sensitive to the phase transition [64].