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

Authors: A. Tawfik
Journal: Advances in High Energy Physics
 (2013)

---

## Abstract

We calculate the first six nonnormalized 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, so that it is possible to evaluate them at chemical freeze-out curve. It is found that a nonmonotonic behaviour 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 the physical observables, we are able to find out the location of the QCD critical endpoint atμ~350 MeV and temperature T~162 MeV.


---

## Body

## 1. 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 cross-over boundary at low density [2, 3]. It is clear that the hadron resonance gas partition function is an approximation to a nonsingular 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–6]. The characteristic feature of CEP is critical dynamical fluctuations [7–14]. The higher order moments are conjectured to reflect the large fluctuations associating the hadron-quark phase transition. This was the motivation of a remarkable number of experimental and theoretical studies [1, 15–20]. 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 ξ [21–23], 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 ofTc≃150-190MeV. As given in [3, 24, 25], depending on the different parameters (for instance, the quark favours and their masses) and on the lattice configurations, lattice QCD assigns different values to Tc. Furthermore, the value of Tc depends on the baryon chemical potential μ. With vanishing μ, the lattice QCD calculations [26] show that the higher order susceptibilities can be related to the higher order moments of the corresponding multiplicity distributions. They also show nonmonotonic behavior near Tc. Apparently, the QCD phase diagram including CEP can be mapped out using this characteristic behaviour. In the literature, one finds that the QCD phase diagram has been determined using various criteria, such as critical energy density [27–29] 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 [30] 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, 31, 32]. 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 Tc [2]. The possibility of finding QCD CEP has motivated the RHIC beam energy scan program [33]. By tuning the collision energy from a center-of-mass energy of 200GeV 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–6, 34]. The higher order moments of various physical quantities are conjectured to reflect the nonmonotonic 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, 27–29, 35, 36] and has been successfully utilized to characterize the conditions deriving from the chemical freeze-out at finite densities [37, 38]. 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 [39]; 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, 27–29, 35, 36]. Such a way, the singularity expected at the Hagedorn temperature [27, 28] 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, Tc.In this paper, we study the first six nonnormalized higher order moments of the particle multiplicity in the hadron resonance gas (HRG) model, which is introduced in Section2. In Section 2.1, 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 2.2. Therefore, it is possible to evaluate them at the chemical freeze-out curve, which is characterized by constant s/T3 at all values of μ, where s and T are entropy density and temperature, respectively. Sections 3 and 4 are devoted to introduce novel conditions describing the chemical freeze-out curve and define the location of the QCD critical endpoint.


## 2. Higher Order Moments of the Particle Multiplicity

In a grand canonical ensemble, it is straightforward to derive an expression for the pressure. In Section2.1, 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 2.2, 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 2.3.The hadron resonances treated as a free gas [3, 27–29, 35] are conjectured to add to the thermodynamic pressure in the hadronic phase (below Tc). 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, 40]. Therefore, the confined phase of QCD, the hadronic phase, is modelled as a noninteracting gas of resonances. The grand canonical partition function reads
(1)Z(T,μ,V)=Tr[exp(μN-H)/T],
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 of the partition function can be calculated exactly and be expressed as a sum over single-particle partition functions Zi1 of all hadrons and their resonances. Consider the following:
(2)lnZ(T,μi,V)=∑i‍lnZi1(T,V)=∑i‍±Vgi2π2∫0∞‍k2dkln{1±exp[μi-εiT]},
where ϵi(k)=(k2+mi2)1/2 is the ith particle dispersion relation, gi 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 [41]. Based on statistical models like Hagedorn [42, 43] and Bootstrap [44, 45], 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 Tc 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 Tc) changed this physical picture drastically. Instead of releasing additional degrees of freedom at T>Tc, the interacting system reduces its effective degrees of freedom at T<Tc. In other words, the hadron gas has much more degrees of freedom than QGP.At finite temperatureT and baryon chemical potential μi, the pressure of the ith hadron or resonance species reads
(3)p(T,μi)=±gi2π2T∫0∞‍k2dkln{1±exp[μi-εiT]}.
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 thehadronic chemical potentials and the quark constituents; μi=3nbμq+nsμS, where nb(ns) 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].


### 2.1. Nonnormalized 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 [46]. Recently, the higher order moments of various multiplicity distributions have been reported by STAR measurements [15–19] and lattice QCD calculations [47–49]. The empirical relevance of the higher order moments to the experimental measurements has been suggested in [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 nonnormalized and normalized moments. The motivation of using resonance gas is clear [27–29]. As introduced about five decades ago, the resonance gas is only valid below Tc [3, 27–29, 35, 42, 43]. Therefore, the present work is valid in the hadronic phase, only.For theith particle, the “first” order moment is given by the derivative of p=-T∂lnZ(T,V,μi)/∂V with respect to the dimensionless quantity μi/T. When taking into account the antiparticles, we add a negative sign to the chemical potential. The first derivative describes the multiplicity distribution or an expectation operator, which is utilized to estimate the number or multiplicity density
(4)m1(T,μi)=±gi2π2T∫0∞‍e(μi-εi)/Tk2dk1±e(μi-εi)/T.
The “second” order moment is known as the variance. It gives the susceptibility of the measurements. Consider the following:
(5)m2(T,μi)=±gi2π2T∫0∞‍e(μi-εi)/Tk2dk1±e(μi-εi)/Tm2(T,μi)=-gi2π2T∫0∞‍e2(μi-εi)/Tk2dk(1±e(μi-εi)/T)2.
The “third” order moment measures the lopsidedness of the distribution. As given in Section 2.2, the normalization of third order moment is known as skewness or the asymmetry in the distribution. Skewness tells us about the direction of variation of the data set. Consider the following:
(6)m3(T,μi)=±gi2π2T∫0∞‍e(μi-εi)/Tk2dk1±e(μi-εi)/Tm3(T,μi)=-gi2π23T∫0∞‍e2(μi-εi)/Tk2dk(1±e(μi-εi)/T)2m3(T,μi)=±gi2π22T∫0∞‍e3(μi-εi)/Tk2dk(1±e(μi-εi)/T)3.
In general, the normalization of rth order moment is obtained by dividing it by σr, where σ is the standard deviation. The normalization is assumed to remove the brightness dependence.The “fourth” order moment compares the tallness and skinniness or shortness and squatness, that is, shape, of a certain measurement to its normal distribution. It is defined as the uncertainty is an uncertainty or“the location- and scale-free movement of probability mass from the shoulders of a distribution into its center and tails and to recognize that it can be formalized in many ways” [50]. Consider the following:
(7)m4(T,μi)=±gi2π2T∫0∞‍e(μi-εi)/Tk2dk1±e(μi-εi)/Tm4(T,μi)=-gi2π27T∫0∞‍e2(μi-εi)/Tk2dk(1±e(μi-εi)/T)2m4(T,μi)=±gi2π212T∫0∞‍e3(μi-εi)/Tk2dk(1±e(μi-εi)/T)3m4(T,μi)=-gi2π26T∫0∞‍e4(μi-εi)/Tk2dk(1±e(μi-εi)/T)4.
The “fifth” order moment measures the asymmetry sensitivity of the “fourth” order moment. Consider the following:
(8)m5(T,μi)=±gi2π2T∫0∞‍e(μi-εi)/Tk2dk1±e(μi-εi)/Tm5(T,μi)=-gi2π215T∫0∞‍e2(μi-εi)/Tk2dk(1±e(μi-εi)/T)2m5(T,μi)=±gi2π250T∫0∞‍e3(μi-εi)/Tk2dk(1±e(μi-εi)/T)3m5(T,μi)=-gi2π260T∫0∞‍e4(μi-εi)/Tk2dk(1±e(μi-εi)/T)4m5(T,μi)=±gi2π224T∫0∞‍e5(μi-εi)/Tk2dk(1±e(μi-εi)/T)5.
The “sixth” order moment is generally associated with compound options. Consider the following:
(9)m6(T,μi)=±gi2π2T∫0∞‍e(μi-εi)/Tk2dk1±e(μi-εi)/Tm6(T,μi)=-gi2π231T∫0∞‍e2(μi-εi)/Tk2dk(1±e(μi-εi)/T)2m6(T,μi)=±gi2π2180T∫0∞‍e3(μi-εi)/Tk2dk(1±e(μi-εi)/T)3m6(T,μi)=-gi2π2390T∫0∞‍e4(μi-εi)/Tk2dk(1±e(μi-εi)/T)4m6(T,μi)=±gi2π2360T∫0∞‍e5(μi-εi)/Tk2dk(1±e(μi-εi)/T)5m6(T,μi)=-gi2π2120T∫0∞‍e6(μi-εi)/Tk2dk(1±e(μi-εi)/T)6.Thus, from (4)–(9), a general expression for the rth order moment can be deduced
(10)mr(T,μi)=gi2π2T∑l=1r‍ar,l∫0∞‍el(μi-εi)/T(1±e(μi-εi)/T)lk2dk,
where the coefficients read
(11)ar,l=(±1)l(-1)l+1[lar-1,l+(l-1)ar-1,l-1],
where l≤r and ar,l vanishes for all 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 [46], 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 rth order. According to [46] and when neglecting the two-particle correlations C2, then the higher order moments read
(12)m2=〈(δN)2〉≈2〈N〉,m3=〈(δN)3〉≈4〈N〉+C3,(13)m4≃〈(δN)4〉=6m3+3m22-8m2+8C4,
where δN=N-〈N〉,
(14)C3(p1,p2,p3)=∑p1‍〈p1p1p1〉+3∑p1<p2‍〈p1p2p2〉C3(p1,p2,p3)=+6∑p1<p2<p3‍〈p1p2p3〉,(15)C4(p1,p2,p3,p4)=∑p1<p2‍〈p1p2p2p2〉+∑p1<p2<p3‍〈p1p2p3p3〉C4(p1,p2,p3,p4)=+8∑p1<p2<p3<p4‍〈p1p2p3p4〉,
and pi is the ith particle. To the second order moment m2, we have to add the effects of the two particle correlation function, 2∑‍C2. The third order moment m3 gets approximately three times this amount. The three and four particle correlations, (14) and (15) appear first in third and fourth moments (13) and (14), 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 [51–57]. In addition to this list of literature, the comprehensive review [58, 59] can be recommended.Taking into account the particle multiplicities, then expressions (12) and (13) can be rewritten as follows:
(16)m2=〈(δN)2〉≃〈N2〉-〈N〉2,m3=〈(δN)3〉≃〈N3〉-〈N2〉〈N〉+2〈N3〉3,m4=〈(δN)4〉-3〈(δN)2〉2≃〈〈N〉4〉-2〈〈N2〉2〉m4=-5〈〈N2〉〉2+6〈〈N〉2〉〈〈N2〉〉.
As mentioned in Section 2.1, the experimental signals for the higher order moments shall be discussed in a future work.As given in (10), the dependence of first six nonnormalized moments on lower ones can explicitly be deduced as follows:
(17)m1=gi2π2T∫0∞‍e(μi-εi)/Tk2dk1±e(μi-εi)/T,m2=m1-gi2π2T∫0∞‍e2(μi-εi)/Tk2dk(1±e(μi-εi)/T)2,m3=-2m1+3m2+2gi2π2T∫0∞‍e3(μi-εi)/Tk2dk(1±e(μi-εi)/T)3,m4=6m1-11m2+6m3-6gi2π2T∫0∞‍e4(-ϵi+μi)/Tk2dk(1±e(-ϵi+μi)/T)4,m5=-24m1+50m2-35m3+10m4m5=+24gi2π2T∫0∞‍e5(-ϵi+μi)/Tk2dk(1±e(-ϵi+μi)/T)5,m6=120m1-274m2+225m3-85m4+15m5m6=-120gi2π2T∫0∞‍e6(-ϵi+μi)/Tk2dk(1±e(-ϵi+μi)/T)6.
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 [46]. From (17), 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(18)mr=(-)r-1(r-1)cm1r-1m1mr=-[(r-1)cm2r-1+(r-2)!]m2+𝒪(m>2),
where cmr is the coefficient of the rth moment. The last term in (18) can be elaborated when higher order moments are calculated. The latter are essential in order to find out a clear pattern.Figure1 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 (18). The dimensionless nonnormalized moments are given depending on the temperature T in the left panel of Figure 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 [46]. 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.(a) Depicts dimensionless nonnormalized moments depending on the temperature. (b) Shows the evolution of the so-called correlation parts (the integrals appearing in (17)) with the temperature. (b) Shows the new coefficients that resulted in when going from lower to higher order moments of multiplicities in a log scale.
(a)
(b)In order to draw a picture about the contributions of the integrals appearing in (17), the phase space integral can be replaced by a series representation. For the ith particle species, the pressure can be expressed as
(19)P(T,μi)=gi2π2T4∑n=1∞‍(±)n+1en(μi/T)(nmiT)2P(T,μi)=×n-4K2(nmiT),
where K2 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 (4), for instance. The resulting expression is applicable as long as μi<mi is valid. Then, the derivatives with respect to μ/T leads to
(20)mr(T,μi)=gi2π2T4∑n=1∞‍(±)n+1en(μi/T)(nmiT)2mr(T,μi)=×nr-4K2(nmiT),
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 the literature. In this case, the expression (19) can be rewritten as
(21)P(T,μi)=giπ2T4∑n=1∞‍(±)n+1(nmiT)2n-4P(T,μi)=×K2(nmiT)cosh(nμiT),
and accordingly,
(22)mr(T,μi)=gi2π2T4∑n=1∞‍(±)n+1(nmiT)2n-4mr(T,μi)=×K2(nmiT)∂r∂μrcosh(nμiT).
Comparing (20) with (19), 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.


### 2.2. 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 a tool to relate moments with various orders to the experimental measurement. The susceptibility of the distribution gives 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 Figure 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 freedom conserved in HRG is introduced in [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 Figure 16. Another feature in these calculations is the assumption that the freeze-out boundary is determined by constant s/T3, where s is the entropy density [37]. At chemical freeze-out boundary, the dependence of σ on μ is given in Figure 5.Results forσ in hadronic (a), bosonic (b), and fermionic (c) resonances are given in dependence on T for various baryon chemical potentials (given in MeV).
(a)
(b)
(c)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 Tc. The QCD CEP is conjectured to be sensitive to skewness. Experimentally, it has been shown that the skewness S is related to ~ξ4.5 [21]. The skewness for bosonic and fermionic resonance gas, respectively, reads
(23)Sb=-12πgiT3/2Sb=×∫0∞‍csch[(εi-μi)/2T]4sinh[(εi-μi)/T]k2dk[∫0∞‍(cosh[(εi-μi)/T]-1)-1k2dk]3/2,Sf=8πgiT3/2Sf=×∫0∞‍csch[(εi-μi)/T]3sinh[(εi-μi)/(2T)]4k2dk[∫0∞‍(cosh[(εi-μi)/T]+1)-1k2dk]3/2.
At different chemical potentials, the skewness S is calculated depending on T and given in Figure 3. It is worthwhile to notice that the values of S in the fermionic resonance gas are much larger than the S values in the bosonic resonance gas, especially at small chemical potentials. But at large chemical potential, this observation is exactly the opposite. The value of μ, at which the large (or small) skewness is flipped, seems to reflect the nature of phase transition, where the fermionic and/or bosonic degrees of freedom turn out to be undistinguishable. Such an observation shall be utilized in positioning the QCD CEP, Section 3.At various baryon chemical potentials (given in MeV), the skewnessS for hadronic (a), bosonic (b), and fermionic (c) resonance gas is given as function of T.
(a)
(b)
(c)The normalization of fourth order moment is known as, heteroskedacity or kurtosis. It means varying volatility or more accurately, varying variance. Actually, the kurtosis is given by normalized fourth order moment minus3. The subtraction of 3, which arises from the Gaussian distribution, is usually omitted [50, 60, 61]. Therefore, the kurtosis is an ideal quantity for probing the non-Gaussian fluctuation feature as expected near Tc and CEP. A sign change of skewness or kurtosis is conjectured to indicate that the system crosses the phase boundary [26, 62, 63]. As HRG is valid below Tc, the sign change is not accessible. It has been shown that kurtosis κ is related to ~ξ7 [21]. Consider the following:
(24)κb=-π2giT3×(([∫0∞‍(1-Cosh[εi-μiT])-1k2dk]2)-1∫0∞‍{cosh[εi-μiT]+2}csch[εi-μi2T]4k2dk×([∫0∞‍(1-cosh[εi-μiT])-1k2dk]2)-1)-3,κf=π2giT3×(([∫0∞‍(cosh[εi-μiT]+1)-1k2dk]2)-1∫0∞‍{cosh[εi-μiT]-2}sech[εi-μi2T]4k2dkmmmm×([∫0∞‍(cosh[εi-μiT]+1)-1k2dk]2)-1)-3.


### 2.3. 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, that is, 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 [26, 64, 65]. Seeking for simplicity, we start with the Boltzmann approximation.


#### 2.3.1. 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(25)σ2〈N〉=∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-2K2(n(mi/T))∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-3K2(n(mi/T)),Sσ=∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-1K2(n(mi/T))∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-2K2(n(mi/T)),κσ2=∑n=1∞(±)n+1en(μi/T)(n(mi/T))2K2(n(mi/T))∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-2K2(n(mi/T))κσ2=-3gi2π2T4∑n=1∞‍(±)n+1en(μi/T)κσ2=×(nmiT)2n-2K2(nmiT),κσS=((∑n=1∞‍(±)n+1en(μi/T)(nmiT)2n-2K2(nmiT))222π2∑n=1∞‍(±)n+1en(μi/T)(nmiT)2K2(nmiT)-3giT4×(∑n=1∞‍(±)n+1en(μi/T)(nmiT)2n-2K2(nmiT))2)×(2π2∑n=1∞‍(±)n+1en(μi/T)(nmiT)2n-1K2(nmiT))-1.
The origin of the second term in all expressions containing κ is obvious. Expressions (25) justify the conclusions in [66] that are in Boltzmann approximation
(26)σ2〈N〉=Sσ≃1,
while
(27)κσ2=κσS≃1-3gi2π2T4exp[μiT](miT)2K2(miT).
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 2.3.2) shows that the products of moments result in a constant dependence on T [66].


#### 2.3.2. 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. The most proposed fluctuation of observables is variations of second order moments of the distribution, such as particle ratio [11–14] and charged dynamical measurement [67]. Then, the fluctuations are approximately related to ξ2 [23].The ratio of standard deviationσ2 and the mean multiplicity 〈N〉 for fermions and bosons read
(28)σ2〈N〉=12∫0∞‍(1±csch[(εi-μi)/T])-1k2dk∫0∞‍(1±e(εi-μi)/T)-1k2dk,
where ± stands for fermions and bosons, respectively. The results are given in Figure 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/T3=7. The ratio σ2/〈N〉 is equivalent to m2/m1, (32). The dependence of bosonic and fermionic σ on μ is given in Figure 5, as well. We notice that σ is smaller than σ2/〈N〉, especially at small μ. At large μ both quantities are almost equal.The ratioσ2/〈N〉 is given depending on T at different chemical potentials μ for hadronic (a), bosonic (b), and fermionic (c) resonance gas.
(a)
(b)
(c)Figure 5
The ratioσ2/〈N〉 and σ are given depending on μ for boson and fermion resonances.The multiplication of skewnessS 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
(29)(Sσ)b=-14∫0∞‍csch[(εi-μi)/T]4sinh[(εi-μi)/T]k2dk∫0∞‍(1-cosh[(εi-μi)/T])-1k2dk,(Sσ)f=4∫0∞‍csch[(εi-μi)/T]3sinh[(εi-μi)/(2T)]4k2dk∫0∞‍(cosh[(εi-μi)/T]+1)-1k2dk.
The product Sσ is equivalent to m3/m2 (33). The results are given in Figure 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 productSσ is given depending on T at different μb values (given in MeV).
(a)
(b)
(c)The dependence ofσ and σ2/〈N〉 on T is illustrated in Figures 2 and 4, respectively. It is obvious that both quantities have a monotonic behavior. Their dependences on μ are given in Figure 5. Also, this type of dependences seems to be monotonic. As given in Figure 6, the product Sσ has a characteristic dependence on T. Regardless of 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 3. When S is studied as a function of μ and given in Figure 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〉, Figure 5. So far, we conclude that starting from the third normalized moment, a nonmonotonic behavior appears. Thus, the kurtosis and its products should be expected to highlight such a nonmonotonic behavior [26, 64, 68].The multiplication of kurtosis byσ2 called κeff [69] is apparently equivalent to the ratio of third order moment to second order moment (34). In lattice QCD and QCD-like models, κeff is found to diverge near the critical endpoint [26, 64]. In HRG, the bosonic and fermionic products read
(30)(κσ2)b=-14((∫0∞‍(1-cosh[εi-μiT])-1k2dk)-1∫0∞‍{cosh[εi-μiT]+2}mmmmmmm×csch[εi-μi2T]4k2dkmmmmmmm×(∫0∞‍(1-cosh[εi-μiT])-1k2dk)-1)+3gi4π21T3∫0∞‍(1-cosh[εi-μiT])-1k2dk,(κσ2)f=14((∫0∞‍(cosh[εi-μiT]+1)-1k2dk)-1∫0∞‍{cosh[εi-μiT]-2}mmmmmmm×sech[εi-μi2T]4k2dkmmmmmmm×(∫0∞‍(cosh[εi-μiT]+1)-1k2dk)-1)-3gi4π21T3∫0∞‍(cosh[εi-μiT]+1)-1k2dk.
When ignoring the constant term in (24), then the second terms in the previous expressions entirely disappear. The results are given in Figure 7. In the hadronic sector, the dependence of κσ2 on the temperature T at different μ values is given in Figure 7(a). 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σ, Figure 6, we simply find that the latter is much drastically changing than the former. Also, when comparing their dependence on the chemical potentials at the freeze-out boundary, Figures 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 4.The productκσ2 is given as a function of T at various μ-values for hadronic (a), bosonic (b), and fermionic (c) resonance gas.
(a)
(b)
(c)The skewness and kurtosis can be combined through the standard deviation as follows:(31)(κσS)b=((∫0∞‍(1-cosh[εi-μiT])-1k2dk)2-14π2T3∫0∞‍{2+cosh[εi-μiT]}mmmmmmm×csch[εi-μi2T]4k2dkmmmmmmm+34gi(∫0∞‍(1-cosh[εi-μiT])-1k2dk)2)(κσS)b=×(14π2T3∫0∞‍csch[εi-μi2T]4mmmmmmm×sinh[εi-μiT]k2dk14π2T3∫0∞‍csch[εi-μi2T]4)-1,(κσS)f=((∫0∞‍(1+cosh[εi-μiT])-1k2dk)214π2T3∫0∞‍{cosh[εi-μiT]-2}mmmmmmn×sech[εi-μi2T]4k2dkmmmmmmn-34gi(∫0∞‍(1+cosh[εi-μiT])-1k2dk)2)(κσS)f=×(4π2T3∫0∞‍csch[εi-μiT]3mmmmmmm×sinh[εi-μi2T]4k2dk∫0∞‍csch[εi-μiT]3)-1.
Also when ignoring the constant term in (24), the second terms in previous expressions disappear. The results of κσ/S are given in Figure 8. In the fermionic sector, κσ/S decreases with increasing T. Increasing μ makes the decrease much faster. Despite the drastic change in κσ/S, this behavior can be compared with Sσ, Figure 6, qualitatively. In the bosonic sector, there is a very slow decrease with increasing μ. The same sector in Sσ, Figure 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 μ. Figure 9 illustrates the μ dependence, that is, κσ/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.κ
σ
/
S as a function T at different chemical potentials for hadrons (a), bosons (b), and fermions (c).
(a)
(b)
(c)Figure 9
The dependence ofκσ/S on the chemical potential μ.From (28)–(31), we summarize that
(32)σ2〈N〉≡m2m1,(33)Sσ≡m3m2,(34)κσ2≡m4m2,(35)κσS≡m4m3.
Obviously, other products would complete missing ratios, for example,
(36)Sχ≡m3m1,κχ2〈N〉≡m4m1,
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 ofSχ≡Sσ2 are given in Figure 10. At the freeze-out boundary, the μ dependence of Sσ2 seems to be much stronger than the μ dependence on Sσ, Figure 11.The thermal evolution ofSσ2 for hadrons (a), bosons (b), and fermions (c) at different μ-values (given in MeV).
(a)
(b)
(c)The dependence ofSσ (a) and Sσ2 (b) of bosons and fermions on the chemical potential μ (given in MeV).
(a)
(b)The thermal evolution ofκσ4/〈N〉 is illustrated in Figure 12. Increasing T is accompanied by an almost overall increase of κσ4/〈N〉 values. At low T, fluctuations appear. In bosonic sector, increasing μ leads to an obvious decrease in κσ4/〈N〉 so that the curves are distributed over a wide range. This is not the case for the fermionic sector. Figure 12(a) shows the results in the hadronic sector. At the freeze-out boundary, the values of κσ4/〈N〉 at various μ-values are given in Figure 13. At large μ, both bosonic and fermionic ratios vanish. With decreasing μ, there is a bounce. It reaches its maximum over ~0.25 GeV. Then both curves come closer to each other. They cross again over ~0.2 GeV. Continue decreasing μ results in another bounce off. At smaller μ, the fermionic sector switches again to positive values, while the bosonic sector results in large negative values of κσ4/〈N〉. For hadrons, κσ4/〈N〉 remains almost negative. At very low μ, it turns out to small positive values.Thermal evolution ofκσ4/〈N〉 for hadrons (a), bosons (b) and fermions (c) at different chemical potentials (given in MeV).
(a)
(b)
(c)Figure 13
The ratioκσ4/〈N〉 is given in depending on μ for boson and fermion resonances.(a) illustrates the chemical freeze-out boundary. The experimental data are given by the solid circles. The curves represent the results of hadron, boson, and fermion gas determined when the sign ofκσ2 is flipped. (b) gives the corresponding s/T3 values. The lattice QCD endpoint is marked by solid square.
(a)
(b)Left panel: Skewness and kurtosis at the freeze-out boundary for fermions and bosons. The curves of both quantities are crossing at almost the sameμ-value. The corresponding T-values are almost equal.
(a)
(b)Figure 16
μ
S as a function of μb in hadron resonance gas model (dotted curve). The linear fitting is given by solid curve.


## 2.1. Nonnormalized 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 [46]. Recently, the higher order moments of various multiplicity distributions have been reported by STAR measurements [15–19] and lattice QCD calculations [47–49]. The empirical relevance of the higher order moments to the experimental measurements has been suggested in [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 nonnormalized and normalized moments. The motivation of using resonance gas is clear [27–29]. As introduced about five decades ago, the resonance gas is only valid below Tc [3, 27–29, 35, 42, 43]. Therefore, the present work is valid in the hadronic phase, only.For theith particle, the “first” order moment is given by the derivative of p=-T∂lnZ(T,V,μi)/∂V with respect to the dimensionless quantity μi/T. When taking into account the antiparticles, we add a negative sign to the chemical potential. The first derivative describes the multiplicity distribution or an expectation operator, which is utilized to estimate the number or multiplicity density
(4)m1(T,μi)=±gi2π2T∫0∞‍e(μi-εi)/Tk2dk1±e(μi-εi)/T.
The “second” order moment is known as the variance. It gives the susceptibility of the measurements. Consider the following:
(5)m2(T,μi)=±gi2π2T∫0∞‍e(μi-εi)/Tk2dk1±e(μi-εi)/Tm2(T,μi)=-gi2π2T∫0∞‍e2(μi-εi)/Tk2dk(1±e(μi-εi)/T)2.
The “third” order moment measures the lopsidedness of the distribution. As given in Section 2.2, the normalization of third order moment is known as skewness or the asymmetry in the distribution. Skewness tells us about the direction of variation of the data set. Consider the following:
(6)m3(T,μi)=±gi2π2T∫0∞‍e(μi-εi)/Tk2dk1±e(μi-εi)/Tm3(T,μi)=-gi2π23T∫0∞‍e2(μi-εi)/Tk2dk(1±e(μi-εi)/T)2m3(T,μi)=±gi2π22T∫0∞‍e3(μi-εi)/Tk2dk(1±e(μi-εi)/T)3.
In general, the normalization of rth order moment is obtained by dividing it by σr, where σ is the standard deviation. The normalization is assumed to remove the brightness dependence.The “fourth” order moment compares the tallness and skinniness or shortness and squatness, that is, shape, of a certain measurement to its normal distribution. It is defined as the uncertainty is an uncertainty or“the location- and scale-free movement of probability mass from the shoulders of a distribution into its center and tails and to recognize that it can be formalized in many ways” [50]. Consider the following:
(7)m4(T,μi)=±gi2π2T∫0∞‍e(μi-εi)/Tk2dk1±e(μi-εi)/Tm4(T,μi)=-gi2π27T∫0∞‍e2(μi-εi)/Tk2dk(1±e(μi-εi)/T)2m4(T,μi)=±gi2π212T∫0∞‍e3(μi-εi)/Tk2dk(1±e(μi-εi)/T)3m4(T,μi)=-gi2π26T∫0∞‍e4(μi-εi)/Tk2dk(1±e(μi-εi)/T)4.
The “fifth” order moment measures the asymmetry sensitivity of the “fourth” order moment. Consider the following:
(8)m5(T,μi)=±gi2π2T∫0∞‍e(μi-εi)/Tk2dk1±e(μi-εi)/Tm5(T,μi)=-gi2π215T∫0∞‍e2(μi-εi)/Tk2dk(1±e(μi-εi)/T)2m5(T,μi)=±gi2π250T∫0∞‍e3(μi-εi)/Tk2dk(1±e(μi-εi)/T)3m5(T,μi)=-gi2π260T∫0∞‍e4(μi-εi)/Tk2dk(1±e(μi-εi)/T)4m5(T,μi)=±gi2π224T∫0∞‍e5(μi-εi)/Tk2dk(1±e(μi-εi)/T)5.
The “sixth” order moment is generally associated with compound options. Consider the following:
(9)m6(T,μi)=±gi2π2T∫0∞‍e(μi-εi)/Tk2dk1±e(μi-εi)/Tm6(T,μi)=-gi2π231T∫0∞‍e2(μi-εi)/Tk2dk(1±e(μi-εi)/T)2m6(T,μi)=±gi2π2180T∫0∞‍e3(μi-εi)/Tk2dk(1±e(μi-εi)/T)3m6(T,μi)=-gi2π2390T∫0∞‍e4(μi-εi)/Tk2dk(1±e(μi-εi)/T)4m6(T,μi)=±gi2π2360T∫0∞‍e5(μi-εi)/Tk2dk(1±e(μi-εi)/T)5m6(T,μi)=-gi2π2120T∫0∞‍e6(μi-εi)/Tk2dk(1±e(μi-εi)/T)6.Thus, from (4)–(9), a general expression for the rth order moment can be deduced
(10)mr(T,μi)=gi2π2T∑l=1r‍ar,l∫0∞‍el(μi-εi)/T(1±e(μi-εi)/T)lk2dk,
where the coefficients read
(11)ar,l=(±1)l(-1)l+1[lar-1,l+(l-1)ar-1,l-1],
where l≤r and ar,l vanishes for all 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 [46], 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 rth order. According to [46] and when neglecting the two-particle correlations C2, then the higher order moments read
(12)m2=〈(δN)2〉≈2〈N〉,m3=〈(δN)3〉≈4〈N〉+C3,(13)m4≃〈(δN)4〉=6m3+3m22-8m2+8C4,
where δN=N-〈N〉,
(14)C3(p1,p2,p3)=∑p1‍〈p1p1p1〉+3∑p1<p2‍〈p1p2p2〉C3(p1,p2,p3)=+6∑p1<p2<p3‍〈p1p2p3〉,(15)C4(p1,p2,p3,p4)=∑p1<p2‍〈p1p2p2p2〉+∑p1<p2<p3‍〈p1p2p3p3〉C4(p1,p2,p3,p4)=+8∑p1<p2<p3<p4‍〈p1p2p3p4〉,
and pi is the ith particle. To the second order moment m2, we have to add the effects of the two particle correlation function, 2∑‍C2. The third order moment m3 gets approximately three times this amount. The three and four particle correlations, (14) and (15) appear first in third and fourth moments (13) and (14), 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 [51–57]. In addition to this list of literature, the comprehensive review [58, 59] can be recommended.Taking into account the particle multiplicities, then expressions (12) and (13) can be rewritten as follows:
(16)m2=〈(δN)2〉≃〈N2〉-〈N〉2,m3=〈(δN)3〉≃〈N3〉-〈N2〉〈N〉+2〈N3〉3,m4=〈(δN)4〉-3〈(δN)2〉2≃〈〈N〉4〉-2〈〈N2〉2〉m4=-5〈〈N2〉〉2+6〈〈N〉2〉〈〈N2〉〉.
As mentioned in Section 2.1, the experimental signals for the higher order moments shall be discussed in a future work.As given in (10), the dependence of first six nonnormalized moments on lower ones can explicitly be deduced as follows:
(17)m1=gi2π2T∫0∞‍e(μi-εi)/Tk2dk1±e(μi-εi)/T,m2=m1-gi2π2T∫0∞‍e2(μi-εi)/Tk2dk(1±e(μi-εi)/T)2,m3=-2m1+3m2+2gi2π2T∫0∞‍e3(μi-εi)/Tk2dk(1±e(μi-εi)/T)3,m4=6m1-11m2+6m3-6gi2π2T∫0∞‍e4(-ϵi+μi)/Tk2dk(1±e(-ϵi+μi)/T)4,m5=-24m1+50m2-35m3+10m4m5=+24gi2π2T∫0∞‍e5(-ϵi+μi)/Tk2dk(1±e(-ϵi+μi)/T)5,m6=120m1-274m2+225m3-85m4+15m5m6=-120gi2π2T∫0∞‍e6(-ϵi+μi)/Tk2dk(1±e(-ϵi+μi)/T)6.
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 [46]. From (17), 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(18)mr=(-)r-1(r-1)cm1r-1m1mr=-[(r-1)cm2r-1+(r-2)!]m2+𝒪(m>2),
where cmr is the coefficient of the rth moment. The last term in (18) can be elaborated when higher order moments are calculated. The latter are essential in order to find out a clear pattern.Figure1 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 (18). The dimensionless nonnormalized moments are given depending on the temperature T in the left panel of Figure 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 [46]. 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.(a) Depicts dimensionless nonnormalized moments depending on the temperature. (b) Shows the evolution of the so-called correlation parts (the integrals appearing in (17)) with the temperature. (b) Shows the new coefficients that resulted in when going from lower to higher order moments of multiplicities in a log scale.
(a)
(b)In order to draw a picture about the contributions of the integrals appearing in (17), the phase space integral can be replaced by a series representation. For the ith particle species, the pressure can be expressed as
(19)P(T,μi)=gi2π2T4∑n=1∞‍(±)n+1en(μi/T)(nmiT)2P(T,μi)=×n-4K2(nmiT),
where K2 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 (4), for instance. The resulting expression is applicable as long as μi<mi is valid. Then, the derivatives with respect to μ/T leads to
(20)mr(T,μi)=gi2π2T4∑n=1∞‍(±)n+1en(μi/T)(nmiT)2mr(T,μi)=×nr-4K2(nmiT),
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 the literature. In this case, the expression (19) can be rewritten as
(21)P(T,μi)=giπ2T4∑n=1∞‍(±)n+1(nmiT)2n-4P(T,μi)=×K2(nmiT)cosh(nμiT),
and accordingly,
(22)mr(T,μi)=gi2π2T4∑n=1∞‍(±)n+1(nmiT)2n-4mr(T,μi)=×K2(nmiT)∂r∂μrcosh(nμiT).
Comparing (20) with (19), 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.


## 2.2. 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 a tool to relate moments with various orders to the experimental measurement. The susceptibility of the distribution gives 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 Figure 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 freedom conserved in HRG is introduced in [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 Figure 16. Another feature in these calculations is the assumption that the freeze-out boundary is determined by constant s/T3, where s is the entropy density [37]. At chemical freeze-out boundary, the dependence of σ on μ is given in Figure 5.Results forσ in hadronic (a), bosonic (b), and fermionic (c) resonances are given in dependence on T for various baryon chemical potentials (given in MeV).
(a)
(b)
(c)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 Tc. The QCD CEP is conjectured to be sensitive to skewness. Experimentally, it has been shown that the skewness S is related to ~ξ4.5 [21]. The skewness for bosonic and fermionic resonance gas, respectively, reads
(23)Sb=-12πgiT3/2Sb=×∫0∞‍csch[(εi-μi)/2T]4sinh[(εi-μi)/T]k2dk[∫0∞‍(cosh[(εi-μi)/T]-1)-1k2dk]3/2,Sf=8πgiT3/2Sf=×∫0∞‍csch[(εi-μi)/T]3sinh[(εi-μi)/(2T)]4k2dk[∫0∞‍(cosh[(εi-μi)/T]+1)-1k2dk]3/2.
At different chemical potentials, the skewness S is calculated depending on T and given in Figure 3. It is worthwhile to notice that the values of S in the fermionic resonance gas are much larger than the S values in the bosonic resonance gas, especially at small chemical potentials. But at large chemical potential, this observation is exactly the opposite. The value of μ, at which the large (or small) skewness is flipped, seems to reflect the nature of phase transition, where the fermionic and/or bosonic degrees of freedom turn out to be undistinguishable. Such an observation shall be utilized in positioning the QCD CEP, Section 3.At various baryon chemical potentials (given in MeV), the skewnessS for hadronic (a), bosonic (b), and fermionic (c) resonance gas is given as function of T.
(a)
(b)
(c)The normalization of fourth order moment is known as, heteroskedacity or kurtosis. It means varying volatility or more accurately, varying variance. Actually, the kurtosis is given by normalized fourth order moment minus3. The subtraction of 3, which arises from the Gaussian distribution, is usually omitted [50, 60, 61]. Therefore, the kurtosis is an ideal quantity for probing the non-Gaussian fluctuation feature as expected near Tc and CEP. A sign change of skewness or kurtosis is conjectured to indicate that the system crosses the phase boundary [26, 62, 63]. As HRG is valid below Tc, the sign change is not accessible. It has been shown that kurtosis κ is related to ~ξ7 [21]. Consider the following:
(24)κb=-π2giT3×(([∫0∞‍(1-Cosh[εi-μiT])-1k2dk]2)-1∫0∞‍{cosh[εi-μiT]+2}csch[εi-μi2T]4k2dk×([∫0∞‍(1-cosh[εi-μiT])-1k2dk]2)-1)-3,κf=π2giT3×(([∫0∞‍(cosh[εi-μiT]+1)-1k2dk]2)-1∫0∞‍{cosh[εi-μiT]-2}sech[εi-μi2T]4k2dkmmmm×([∫0∞‍(cosh[εi-μiT]+1)-1k2dk]2)-1)-3.


## 2.3. 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, that is, 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 [26, 64, 65]. Seeking for simplicity, we start with the Boltzmann approximation.


### 2.3.1. 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(25)σ2〈N〉=∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-2K2(n(mi/T))∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-3K2(n(mi/T)),Sσ=∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-1K2(n(mi/T))∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-2K2(n(mi/T)),κσ2=∑n=1∞(±)n+1en(μi/T)(n(mi/T))2K2(n(mi/T))∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-2K2(n(mi/T))κσ2=-3gi2π2T4∑n=1∞‍(±)n+1en(μi/T)κσ2=×(nmiT)2n-2K2(nmiT),κσS=((∑n=1∞‍(±)n+1en(μi/T)(nmiT)2n-2K2(nmiT))222π2∑n=1∞‍(±)n+1en(μi/T)(nmiT)2K2(nmiT)-3giT4×(∑n=1∞‍(±)n+1en(μi/T)(nmiT)2n-2K2(nmiT))2)×(2π2∑n=1∞‍(±)n+1en(μi/T)(nmiT)2n-1K2(nmiT))-1.
The origin of the second term in all expressions containing κ is obvious. Expressions (25) justify the conclusions in [66] that are in Boltzmann approximation
(26)σ2〈N〉=Sσ≃1,
while
(27)κσ2=κσS≃1-3gi2π2T4exp[μiT](miT)2K2(miT).
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 2.3.2) shows that the products of moments result in a constant dependence on T [66].


### 2.3.2. 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. The most proposed fluctuation of observables is variations of second order moments of the distribution, such as particle ratio [11–14] and charged dynamical measurement [67]. Then, the fluctuations are approximately related to ξ2 [23].The ratio of standard deviationσ2 and the mean multiplicity 〈N〉 for fermions and bosons read
(28)σ2〈N〉=12∫0∞‍(1±csch[(εi-μi)/T])-1k2dk∫0∞‍(1±e(εi-μi)/T)-1k2dk,
where ± stands for fermions and bosons, respectively. The results are given in Figure 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/T3=7. The ratio σ2/〈N〉 is equivalent to m2/m1, (32). The dependence of bosonic and fermionic σ on μ is given in Figure 5, as well. We notice that σ is smaller than σ2/〈N〉, especially at small μ. At large μ both quantities are almost equal.The ratioσ2/〈N〉 is given depending on T at different chemical potentials μ for hadronic (a), bosonic (b), and fermionic (c) resonance gas.
(a)
(b)
(c)Figure 5
The ratioσ2/〈N〉 and σ are given depending on μ for boson and fermion resonances.The multiplication of skewnessS 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
(29)(Sσ)b=-14∫0∞‍csch[(εi-μi)/T]4sinh[(εi-μi)/T]k2dk∫0∞‍(1-cosh[(εi-μi)/T])-1k2dk,(Sσ)f=4∫0∞‍csch[(εi-μi)/T]3sinh[(εi-μi)/(2T)]4k2dk∫0∞‍(cosh[(εi-μi)/T]+1)-1k2dk.
The product Sσ is equivalent to m3/m2 (33). The results are given in Figure 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 productSσ is given depending on T at different μb values (given in MeV).
(a)
(b)
(c)The dependence ofσ and σ2/〈N〉 on T is illustrated in Figures 2 and 4, respectively. It is obvious that both quantities have a monotonic behavior. Their dependences on μ are given in Figure 5. Also, this type of dependences seems to be monotonic. As given in Figure 6, the product Sσ has a characteristic dependence on T. Regardless of 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 3. When S is studied as a function of μ and given in Figure 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〉, Figure 5. So far, we conclude that starting from the third normalized moment, a nonmonotonic behavior appears. Thus, the kurtosis and its products should be expected to highlight such a nonmonotonic behavior [26, 64, 68].The multiplication of kurtosis byσ2 called κeff [69] is apparently equivalent to the ratio of third order moment to second order moment (34). In lattice QCD and QCD-like models, κeff is found to diverge near the critical endpoint [26, 64]. In HRG, the bosonic and fermionic products read
(30)(κσ2)b=-14((∫0∞‍(1-cosh[εi-μiT])-1k2dk)-1∫0∞‍{cosh[εi-μiT]+2}mmmmmmm×csch[εi-μi2T]4k2dkmmmmmmm×(∫0∞‍(1-cosh[εi-μiT])-1k2dk)-1)+3gi4π21T3∫0∞‍(1-cosh[εi-μiT])-1k2dk,(κσ2)f=14((∫0∞‍(cosh[εi-μiT]+1)-1k2dk)-1∫0∞‍{cosh[εi-μiT]-2}mmmmmmm×sech[εi-μi2T]4k2dkmmmmmmm×(∫0∞‍(cosh[εi-μiT]+1)-1k2dk)-1)-3gi4π21T3∫0∞‍(cosh[εi-μiT]+1)-1k2dk.
When ignoring the constant term in (24), then the second terms in the previous expressions entirely disappear. The results are given in Figure 7. In the hadronic sector, the dependence of κσ2 on the temperature T at different μ values is given in Figure 7(a). 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σ, Figure 6, we simply find that the latter is much drastically changing than the former. Also, when comparing their dependence on the chemical potentials at the freeze-out boundary, Figures 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 4.The productκσ2 is given as a function of T at various μ-values for hadronic (a), bosonic (b), and fermionic (c) resonance gas.
(a)
(b)
(c)The skewness and kurtosis can be combined through the standard deviation as follows:(31)(κσS)b=((∫0∞‍(1-cosh[εi-μiT])-1k2dk)2-14π2T3∫0∞‍{2+cosh[εi-μiT]}mmmmmmm×csch[εi-μi2T]4k2dkmmmmmmm+34gi(∫0∞‍(1-cosh[εi-μiT])-1k2dk)2)(κσS)b=×(14π2T3∫0∞‍csch[εi-μi2T]4mmmmmmm×sinh[εi-μiT]k2dk14π2T3∫0∞‍csch[εi-μi2T]4)-1,(κσS)f=((∫0∞‍(1+cosh[εi-μiT])-1k2dk)214π2T3∫0∞‍{cosh[εi-μiT]-2}mmmmmmn×sech[εi-μi2T]4k2dkmmmmmmn-34gi(∫0∞‍(1+cosh[εi-μiT])-1k2dk)2)(κσS)f=×(4π2T3∫0∞‍csch[εi-μiT]3mmmmmmm×sinh[εi-μi2T]4k2dk∫0∞‍csch[εi-μiT]3)-1.
Also when ignoring the constant term in (24), the second terms in previous expressions disappear. The results of κσ/S are given in Figure 8. In the fermionic sector, κσ/S decreases with increasing T. Increasing μ makes the decrease much faster. Despite the drastic change in κσ/S, this behavior can be compared with Sσ, Figure 6, qualitatively. In the bosonic sector, there is a very slow decrease with increasing μ. The same sector in Sσ, Figure 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 μ. Figure 9 illustrates the μ dependence, that is, κσ/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.κ
σ
/
S as a function T at different chemical potentials for hadrons (a), bosons (b), and fermions (c).
(a)
(b)
(c)Figure 9
The dependence ofκσ/S on the chemical potential μ.From (28)–(31), we summarize that
(32)σ2〈N〉≡m2m1,(33)Sσ≡m3m2,(34)κσ2≡m4m2,(35)κσS≡m4m3.
Obviously, other products would complete missing ratios, for example,
(36)Sχ≡m3m1,κχ2〈N〉≡m4m1,
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 ofSχ≡Sσ2 are given in Figure 10. At the freeze-out boundary, the μ dependence of Sσ2 seems to be much stronger than the μ dependence on Sσ, Figure 11.The thermal evolution ofSσ2 for hadrons (a), bosons (b), and fermions (c) at different μ-values (given in MeV).
(a)
(b)
(c)The dependence ofSσ (a) and Sσ2 (b) of bosons and fermions on the chemical potential μ (given in MeV).
(a)
(b)The thermal evolution ofκσ4/〈N〉 is illustrated in Figure 12. Increasing T is accompanied by an almost overall increase of κσ4/〈N〉 values. At low T, fluctuations appear. In bosonic sector, increasing μ leads to an obvious decrease in κσ4/〈N〉 so that the curves are distributed over a wide range. This is not the case for the fermionic sector. Figure 12(a) shows the results in the hadronic sector. At the freeze-out boundary, the values of κσ4/〈N〉 at various μ-values are given in Figure 13. At large μ, both bosonic and fermionic ratios vanish. With decreasing μ, there is a bounce. It reaches its maximum over ~0.25 GeV. Then both curves come closer to each other. They cross again over ~0.2 GeV. Continue decreasing μ results in another bounce off. At smaller μ, the fermionic sector switches again to positive values, while the bosonic sector results in large negative values of κσ4/〈N〉. For hadrons, κσ4/〈N〉 remains almost negative. At very low μ, it turns out to small positive values.Thermal evolution ofκσ4/〈N〉 for hadrons (a), bosons (b) and fermions (c) at different chemical potentials (given in MeV).
(a)
(b)
(c)Figure 13
The ratioκσ4/〈N〉 is given in depending on μ for boson and fermion resonances.(a) illustrates the chemical freeze-out boundary. The experimental data are given by the solid circles. The curves represent the results of hadron, boson, and fermion gas determined when the sign ofκσ2 is flipped. (b) gives the corresponding s/T3 values. The lattice QCD endpoint is marked by solid square.
(a)
(b)Left panel: Skewness and kurtosis at the freeze-out boundary for fermions and bosons. The curves of both quantities are crossing at almost the sameμ-value. The corresponding T-values are almost equal.
(a)
(b)Figure 16
μ
S as a function of μb in hadron resonance gas model (dotted curve). The linear fitting is given by solid curve.


## 2.3.1. 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(25)σ2〈N〉=∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-2K2(n(mi/T))∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-3K2(n(mi/T)),Sσ=∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-1K2(n(mi/T))∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-2K2(n(mi/T)),κσ2=∑n=1∞(±)n+1en(μi/T)(n(mi/T))2K2(n(mi/T))∑n=1∞(±)n+1en(μi/T)(n(mi/T))2n-2K2(n(mi/T))κσ2=-3gi2π2T4∑n=1∞‍(±)n+1en(μi/T)κσ2=×(nmiT)2n-2K2(nmiT),κσS=((∑n=1∞‍(±)n+1en(μi/T)(nmiT)2n-2K2(nmiT))222π2∑n=1∞‍(±)n+1en(μi/T)(nmiT)2K2(nmiT)-3giT4×(∑n=1∞‍(±)n+1en(μi/T)(nmiT)2n-2K2(nmiT))2)×(2π2∑n=1∞‍(±)n+1en(μi/T)(nmiT)2n-1K2(nmiT))-1.
The origin of the second term in all expressions containing κ is obvious. Expressions (25) justify the conclusions in [66] that are in Boltzmann approximation
(26)σ2〈N〉=Sσ≃1,
while
(27)κσ2=κσS≃1-3gi2π2T4exp[μiT](miT)2K2(miT).
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 2.3.2) shows that the products of moments result in a constant dependence on T [66].


## 2.3.2. 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. The most proposed fluctuation of observables is variations of second order moments of the distribution, such as particle ratio [11–14] and charged dynamical measurement [67]. Then, the fluctuations are approximately related to ξ2 [23].The ratio of standard deviationσ2 and the mean multiplicity 〈N〉 for fermions and bosons read
(28)σ2〈N〉=12∫0∞‍(1±csch[(εi-μi)/T])-1k2dk∫0∞‍(1±e(εi-μi)/T)-1k2dk,
where ± stands for fermions and bosons, respectively. The results are given in Figure 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/T3=7. The ratio σ2/〈N〉 is equivalent to m2/m1, (32). The dependence of bosonic and fermionic σ on μ is given in Figure 5, as well. We notice that σ is smaller than σ2/〈N〉, especially at small μ. At large μ both quantities are almost equal.The ratioσ2/〈N〉 is given depending on T at different chemical potentials μ for hadronic (a), bosonic (b), and fermionic (c) resonance gas.
(a)
(b)
(c)Figure 5
The ratioσ2/〈N〉 and σ are given depending on μ for boson and fermion resonances.The multiplication of skewnessS 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
(29)(Sσ)b=-14∫0∞‍csch[(εi-μi)/T]4sinh[(εi-μi)/T]k2dk∫0∞‍(1-cosh[(εi-μi)/T])-1k2dk,(Sσ)f=4∫0∞‍csch[(εi-μi)/T]3sinh[(εi-μi)/(2T)]4k2dk∫0∞‍(cosh[(εi-μi)/T]+1)-1k2dk.
The product Sσ is equivalent to m3/m2 (33). The results are given in Figure 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 productSσ is given depending on T at different μb values (given in MeV).
(a)
(b)
(c)The dependence ofσ and σ2/〈N〉 on T is illustrated in Figures 2 and 4, respectively. It is obvious that both quantities have a monotonic behavior. Their dependences on μ are given in Figure 5. Also, this type of dependences seems to be monotonic. As given in Figure 6, the product Sσ has a characteristic dependence on T. Regardless of 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 3. When S is studied as a function of μ and given in Figure 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〉, Figure 5. So far, we conclude that starting from the third normalized moment, a nonmonotonic behavior appears. Thus, the kurtosis and its products should be expected to highlight such a nonmonotonic behavior [26, 64, 68].The multiplication of kurtosis byσ2 called κeff [69] is apparently equivalent to the ratio of third order moment to second order moment (34). In lattice QCD and QCD-like models, κeff is found to diverge near the critical endpoint [26, 64]. In HRG, the bosonic and fermionic products read
(30)(κσ2)b=-14((∫0∞‍(1-cosh[εi-μiT])-1k2dk)-1∫0∞‍{cosh[εi-μiT]+2}mmmmmmm×csch[εi-μi2T]4k2dkmmmmmmm×(∫0∞‍(1-cosh[εi-μiT])-1k2dk)-1)+3gi4π21T3∫0∞‍(1-cosh[εi-μiT])-1k2dk,(κσ2)f=14((∫0∞‍(cosh[εi-μiT]+1)-1k2dk)-1∫0∞‍{cosh[εi-μiT]-2}mmmmmmm×sech[εi-μi2T]4k2dkmmmmmmm×(∫0∞‍(cosh[εi-μiT]+1)-1k2dk)-1)-3gi4π21T3∫0∞‍(cosh[εi-μiT]+1)-1k2dk.
When ignoring the constant term in (24), then the second terms in the previous expressions entirely disappear. The results are given in Figure 7. In the hadronic sector, the dependence of κσ2 on the temperature T at different μ values is given in Figure 7(a). 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σ, Figure 6, we simply find that the latter is much drastically changing than the former. Also, when comparing their dependence on the chemical potentials at the freeze-out boundary, Figures 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 4.The productκσ2 is given as a function of T at various μ-values for hadronic (a), bosonic (b), and fermionic (c) resonance gas.
(a)
(b)
(c)The skewness and kurtosis can be combined through the standard deviation as follows:(31)(κσS)b=((∫0∞‍(1-cosh[εi-μiT])-1k2dk)2-14π2T3∫0∞‍{2+cosh[εi-μiT]}mmmmmmm×csch[εi-μi2T]4k2dkmmmmmmm+34gi(∫0∞‍(1-cosh[εi-μiT])-1k2dk)2)(κσS)b=×(14π2T3∫0∞‍csch[εi-μi2T]4mmmmmmm×sinh[εi-μiT]k2dk14π2T3∫0∞‍csch[εi-μi2T]4)-1,(κσS)f=((∫0∞‍(1+cosh[εi-μiT])-1k2dk)214π2T3∫0∞‍{cosh[εi-μiT]-2}mmmmmmn×sech[εi-μi2T]4k2dkmmmmmmn-34gi(∫0∞‍(1+cosh[εi-μiT])-1k2dk)2)(κσS)f=×(4π2T3∫0∞‍csch[εi-μiT]3mmmmmmm×sinh[εi-μi2T]4k2dk∫0∞‍csch[εi-μiT]3)-1.
Also when ignoring the constant term in (24), the second terms in previous expressions disappear. The results of κσ/S are given in Figure 8. In the fermionic sector, κσ/S decreases with increasing T. Increasing μ makes the decrease much faster. Despite the drastic change in κσ/S, this behavior can be compared with Sσ, Figure 6, qualitatively. In the bosonic sector, there is a very slow decrease with increasing μ. The same sector in Sσ, Figure 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 μ. Figure 9 illustrates the μ dependence, that is, κσ/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.κ
σ
/
S as a function T at different chemical potentials for hadrons (a), bosons (b), and fermions (c).
(a)
(b)
(c)Figure 9
The dependence ofκσ/S on the chemical potential μ.From (28)–(31), we summarize that
(32)σ2〈N〉≡m2m1,(33)Sσ≡m3m2,(34)κσ2≡m4m2,(35)κσS≡m4m3.
Obviously, other products would complete missing ratios, for example,
(36)Sχ≡m3m1,κχ2〈N〉≡m4m1,
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 ofSχ≡Sσ2 are given in Figure 10. At the freeze-out boundary, the μ dependence of Sσ2 seems to be much stronger than the μ dependence on Sσ, Figure 11.The thermal evolution ofSσ2 for hadrons (a), bosons (b), and fermions (c) at different μ-values (given in MeV).
(a)
(b)
(c)The dependence ofSσ (a) and Sσ2 (b) of bosons and fermions on the chemical potential μ (given in MeV).
(a)
(b)The thermal evolution ofκσ4/〈N〉 is illustrated in Figure 12. Increasing T is accompanied by an almost overall increase of κσ4/〈N〉 values. At low T, fluctuations appear. In bosonic sector, increasing μ leads to an obvious decrease in κσ4/〈N〉 so that the curves are distributed over a wide range. This is not the case for the fermionic sector. Figure 12(a) shows the results in the hadronic sector. At the freeze-out boundary, the values of κσ4/〈N〉 at various μ-values are given in Figure 13. At large μ, both bosonic and fermionic ratios vanish. With decreasing μ, there is a bounce. It reaches its maximum over ~0.25 GeV. Then both curves come closer to each other. They cross again over ~0.2 GeV. Continue decreasing μ results in another bounce off. At smaller μ, the fermionic sector switches again to positive values, while the bosonic sector results in large negative values of κσ4/〈N〉. For hadrons, κσ4/〈N〉 remains almost negative. At very low μ, it turns out to small positive values.Thermal evolution ofκσ4/〈N〉 for hadrons (a), bosons (b) and fermions (c) at different chemical potentials (given in MeV).
(a)
(b)
(c)Figure 13
The ratioκσ4/〈N〉 is given in depending on μ for boson and fermion resonances.(a) illustrates the chemical freeze-out boundary. The experimental data are given by the solid circles. The curves represent the results of hadron, boson, and fermion gas determined when the sign ofκσ2 is flipped. (b) gives the corresponding s/T3 values. The lattice QCD endpoint is marked by solid square.
(a)
(b)Left panel: Skewness and kurtosis at the freeze-out boundary for fermions and bosons. The curves of both quantities are crossing at almost the sameμ-value. The corresponding T-values are almost equal.
(a)
(b)Figure 16
μ
S as a function of μb in hadron resonance gas model (dotted curve). The linear fitting is given by solid curve.


## 3. Chemical Freeze-Out

In the rest frame of produced particle, the hadronic matter can be determined by constant degrees of freedom, for instance,s/T3(4/π2)=const [70]. The chemical freeze-out is related to the particle creation. Therefore, the abundance of different particle species is controlled by chemical potential, which obviously depends on T. With the beam energy, T is increasing, while the baryon densities at midrapidity are 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 follow regular patterns as the beam energy increases [30, 37]. The higher order moments have been suggested to control the chemical freeze-out, so that several conditions have been proposed [71–73].As introduced in Section2.3.2, the thermal evolution of κσ2 declines. This result seems to update previous ones [48, 74], 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 [71–73]. It has been found that the T and μ parameters, at which the sign is flipped, amazingly coincide with the ones of the chemical freeze-out, Figure 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 4. The results are given in Figure 14. Vanishing κσ2 leads to
(37)∫0∞‍{cosh[(εi-μi)T]+2}csch[(εi-μi)2T]4k2dk=3giπ21T3[∫0∞‍(1-cosh[(εi-μi)T])-1k2dk]2,∫0∞‍{cosh[(εi-μi)T]-2}sech[(εi-μi)2T]4k2dk=3giπ21T3[∫0∞‍(cosh[(εi-μi)T]+1)-1k2dk]2.
The rhs and lhs in both expressions can be rewritten as
(38)16π2giT3κ=48π2giT3χ2,
which is valid for bosons and fermions. Then, the chemical freeze-out is defined, if the condition
(39)κ(T,μ)=3χ2(T,μ)
is fulfilled. At chemical freeze-out curve, a naive estimation leads to ξ~31/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 [75]. 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 parametersT 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 Figure 14 collects a large experimental data set. For a recent review, we refer to [30] 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 [37]; s/T3=const., right panel of Figure 14, at least over the range 0<μ<0.6 GeV. When excluding all degrees of freedom but the fermionic ones, s/T3 decreases with increasing μ (dashed curve). An opposite dependence is accompanied with the bosonic degrees of freedom (dotted curve).


## 4. QCD Critical Endpoint

In the right panel of Figure14, 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 Figure 14. Furthermore, based on nonperturbative convergence radius, the critical endpoint as calculated in lattice QCD [5, 6] given by solid square in Figure 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 nonsingular 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 Figure 15.In Figure15, the skewness and kurtosis of bosons and fermions calculated at the freeze-out curve are given depending on μ (left panel) and center-of-mass energy s (right panel). It seems that the evolution of fermionic and bosonic skewness and kurtosis coincide at one point, marked with the vertical band. This fact would reflect the nature of the phase transition. It would be the critical endpoint connecting cross-over with the first order deconfinement phase transitions. At QCD CEP, the phase transition is conjectured to be of second order. It is worthwhile to mention that the crossing point is amazingly coincident with the QCD CEP measured in lattice QCD [5, 6], Figure 14, regardless of its uncertainties.From (23) and (24), the following expressions have to be solved in μ, individually and/or dependently, in order to determine μ of the crossing point:
(40)-∫0∞‍csch[(εi-μi)/2T]4sinh[(εi-μi)/T]k2dk[∫0∞‍(cosh[(εi-μi)/T]-1)-1k2dk]3/2=16(([∫0∞‍(cosh[εi-μiT]+1)-1k2dk]3/2)-1∫0∞‍csch[εi-μiT]3sinh[εi-μi2T]4k2dkmmmmmm×([∫0∞‍(cosh[εi-μiT]+1)-1k2dk]3/2)-1),-∫0∞‍{cosh[(εi-μi)/T]+2}csch[(εi-μi)/(2T)]4k2dk[∫0∞‍(1-cosh[(εi-μi)/T])-1k2dk]2=∫0∞‍{cosh[εi-μiT]-2}sech[εi-μi2T]4k2dk×([∫0∞‍(cosh[εi-μiT]+1)-1k2dk]2)-1.
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 2.3.1.
(41)p(T,μq,μs)T4=±12π2T3∑i=1∞‍gimi2∑n=1∞‍(±)n+1n2p(T,μq,μs)T4=×K2(nmiT)exp[n(3nb+ns)μq-nsμsT],
where nb and ns 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 2. 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
(42)μb=3nbμq=Tln[gbmb2K2(mb/T)gfmf2K2(mf/T)].
In the relativistic limit, K2(m/T)≈2T2/m2-1/2 while in the nonrelativistic limit K2(m/T)≈πT/2mexp(-m/T)(1+15T/8m). It is obvious that the bosonic and fermionic degrees of freedom play an essential role in determining (42). 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 Figure 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.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 [76]. 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 [77].


## 5. Conclusion

In the present work, the first six nonnormalized order moments of the particle multiplicity are calculated in the hadron resonance gas model. 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 that consists of all 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/T3 at all values of μ, where s and T are entropy density and temperature, respectively. It has been found that nonmonotonic behavior reflecting dynamical fluctuation and strong correlations appears starting from the normalized third order moment (skewness S). Furthermore, nonmonotonicity is observed in the normalized fourth order moment, the kurtosis κ, and its products. These are novel observations. Although HRG is exclusively applicable below Tc, that is, 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 evidence 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, that is, the second order moment. The location of QCD critical endpoint (at T and μ axes) is positioned when the condition Sb=Sf or κb=κf is fulfilled. The subscripts b and f refer to bosons and fermions, respectively. Accordingly, the QCD endpoint is positioned at μ~350MeV and T~162MeV. 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 was posted in arXiv, in which the authors used the second order moment (susceptibility) and the fourth order one (kurtosis) that are apparently sensitive to the phase transition [78].



---

*Source: 574871-2013-07-18.xml*