3.1 Hadron Resonance Gas Model

In the Hadron Resonance Gas (HRG) model, non-interacting hadrons and their resonance are the basic degree of freedom. The interactions are encoded in the thermal creation of hadronic resonances based on their Boltzmann factor. The HRG can successfully describe the observed particle abundances in heavy ion collisions. For simplify and discussion purpose, we use the Boltzmann approximation and the pressure can be expressed [49-51]

where gi is the degeneracy factor for hadrons of mass mi, and ${\widehat{\mu }}_q\equiv \frac{{\mu }_q}{T}$, with q=B, S, Q denote the net-baryon number, net-strangeness and the net-charge, and μB, μS, μQ are the corresponding chemical potentials respectively. The K₂(x) is the modified Bessel function and the summation is taking over all stable hadrons and resonance and thus the contribution of anti-particles are automatically included. The results from HRG model are usually served as a baseline for finding the signature of phase transition and QCD critical point in heavy-ion collisions. For net-baryon number fluctuations, the ratios of cumulants from HRG model are simple. With Boltzmann approximation, the baryon number susceptibility can be expressed as

Thus, the ratio of baryon number susceptibilities can be easily obtained

Based on the Eq.(15) and (16), we obtain

where n=1,2,3..., μB and T are the baryon chemical potential and temperature of the thermal system. This simple result arises from the fact that only baryons with baryon number B = 1 contribute to the various cumulants in the HRG model. However, due to the contribution of the multi-charge states Q = 2 or S = 2,3 for net-charge and net-strangeness fluctuations, respectively, thus the results of net-charge and net-strangeness fluctuations are more complicated than net-baryon number fluctuations from the HRG model. Figure 4 shows the ratio of susceptibilities of charge (left) and strangeness (right) from the HRG model calculations along the chemical freeze-out curve in heavy-ion collisions. It can be found that the susceptibilities ratios χ₂/χ₁ or χ₃/χ₂ of charge and strangeness show strong energy dependence whereas the χ₄/χ₂ has small variation with energies. Due to the contributions from the multi-charge states Q = 2 or S = 2,3, the charge and strangeness χ₄/χ₂ deviate from unity.

Figure 4:Full-screen View

(Color online) The ratio of susceptibilities of charge (left) and strangeness (right) fluctuations as a function of colliding energy along the parameterized freeze-out curve in heavy-ion collisions. The results are calculated from the HRG model [49].

3.2 Lattice QCD

Lattice QCD is a well-established non-perturbative approach to solve the QCD theory of quarks and gluons exactly from first principles and without any assumptions [52]. It can be used to study the thermodynamic properties of a strongly interacting system in thermal equilibrium. Most importantly, Lattice QCD provides a framework for investigation of non-perturbative phenomena such as confinement and quark-gluon plasma formation, which are intractable by means of analytic field theories.

Figure 5 shows the results of QCD equation of state (the trace anomaly, the pressure and the entropy density) from two independent groups: Hot QCD and Wuppertal-Budapest Collaboration, which used the different actions. The results from the two groups got good agreement with each other. On the other hand, the pressure (P/T⁴) at finite μB region can be calculated by using the Taylor expansion techniques. By putting the μQ=μS=0, we can expand the pressure (P/T⁴) into finite μB as [20, 39, 54]

Figure 5:Full-screen View

(Color online) The comparison of the HISQ/tree (Hot QCD Collaboration) and stout (Wuppertal-Budapest Collaboration) results for the trace anomaly, the pressure, and the entropy density in the Lattice QCD calculation at vanishing baryon chemical potential [53].

Due to the symmetry of QCD, the odd terms are vanishing and only even terms are left. It shows various order corrections to the pressure. The coefficients of leading order (LO), next leading order (NLO) and next next leading order (NNLO) are related to the baryon number susceptibilities ${\chi }_2^B$, ${\chi }_4^B/{\chi }_2^B$ and ${\chi }_6^B/{\chi }_2^B$, respectively. Those susceptibilities are defined in Eq. (2) and can be evaluated at μB=μ_Q=μS=0. Figures 6 and 7 left show the preliminary BNL-Bielefeld-CCNU results of baryon number susceptibility (${\chi }_2^B$) and the susceptibilities ratios (${\chi }_4^B/{\chi }_2^B$ and ${\chi }_6^B/{\chi }_2^B$) as a function of temperature computed from Lattice QCD at vanishing chemical potentials. It can be found that at low temperature, the results from Lattice QCD are consistent with the results from HRG whereas it shows large discrepancies between those two at high temperature. The ratio ${\chi }_6^B/{\chi }_2^B$ shows negative values near and above transition temperature and positive values around unity at low temperatures, but there are still large uncertainties and more statistics are needed. Figure 7 right shows the Lattice calculation of the pressure at finite baryon density and the effects of the correction with different orders. The correction of the NNLO term on the pressure is found to be very small (<5%) when μB/T<2. It means the Taylor expansion up to NNLO order is under control with μB/T<2. The susceptibility of conserved charges (B, Q, S) can be also computed at finite baryon density region by using the Taylor series in terms of the baryon chemical potential (μB) at μQ=μS=0 [20, 39, 54]

Figure 6:Full-screen View

(Color online) Lattice QCD results from BNL-Bielefeld-CCNU collaboration [20, 39, 54]: the second order baryon number susceptibility (${\chi }_2^B$) (left) and the fourth to second order baryon number susceptibilities ratio (${\chi }_4^B/{\chi }_2^B$) (right) as a function of temperature calculated from Lattice QCD at vanishing chemical potential (μq=0, q=B,Q,S)

Figure 7:Full-screen View

(Color online) Lattice QCD results from BNL-Bielefeld-CCNU collaboration [20, 39, 54]: the sixth to second order baryon number susceptibilities ratio (${\chi }_6^B/{\chi }_2^B$) at vanishing chemical potential (left) and the pressure as a function of temperature for different baryon chemical potential and precision level (right).

In the following, we focus on discussing the next leading order (NLO) Taylor expansion of the baryon number susceptibilities. Due to the QCD sysmetry for matter and anti-matter, the NLO Taylor expansion for the odd and even order baryon number susceptibilities can be expressed as

where the ${\chi }_{2n}^B(T)$ and ${\chi }_{2n+2}^B(T)$ are the baryon number susceptibilities evaluated at μB=μQ=μS=0 with n=1,2,3,4,...,N. If we define a dimensionless quantity Ln:

With this definition, the Taylor expansion of the odd and even order susceptibility at the next leading order can be re-written as

Then, we can express the baryon number susceptibilities ratios as

If we consider Ln<<1, the r.h.s. of the Eq. (28) to (31) can be simplified as

Based on the Eqs. (32), (34), and (35), we have

where we use a leading order approximation $\frac{{\sigma }_2^B}{M_B}(T\mathrm{,}{\mu }_B)=\frac{T}{{\mu }_B}$. We define a temperature dependent quantity rn(T):

Then, deriving from Eqs. (36) and (37), we get

For the lowest order with n = 1, we obtain

where $L_1=\frac{1}{6}\frac{{\chi }_4^B(T)}{{\chi }_2^B(T)}{\left(\frac{{\mu }_B}{T}\right)}^2$ and $L_2=\frac{1}{6}\frac{{\chi }_6^B(T)}{{\chi }_4^B(T)}{\left(\frac{{\mu }_B}{T}\right)}^2$. We may find that ${\kappa }_B{\sigma }_B^2$ and $S_B{\sigma }_B^3/M_B$ are closely related and their difference is

where

In above Taylor expansions of the baryon number susceptibilities in Lattice QCD, we always assume μQ=μS=0 and expand up to next to leading order. For more realistic, one need to consider the case μQ ≠μS ≠ 0. In order to compare with experimental data, we need additional constrains. For example, strangeness neutrality (NS=0) and baryon to charge number ratios equals to 0.4 (NQ/NB=0.4) in Au+Au and Pb+Pb collisions. Furthermore, the self-consistent determination of the freeze-out in QCD thermodynamics for heavy-ion collisions are needed and makes the comparison between Lattice QCD and experimental data with more complication [38, 39].

4.1 Fluctuation Signature near QCD Critical Point

The characteristic feature of critical point is the divergence of the correlation length, which is limited by the system size and finite time effects due to the critical slowing down. When the critical point is passed by the thermodynamic condition of the matter created in heavy-ion collisions, the expected signature is the non-monotonic variation of the observables with the colliding energy. Many theoretical and model calculations including critical fluctuations have been done for the fluctuations of conserved charges (B,Q,S) along the chemical freeze-out lines in heavy-ion collisions. Those can provide predictions on the energy dependence of the fluctuation observables when passing by the critical point.

4.1.1 σ Field Model

One of the most important calculations is done with the σ field model [18]. This calculations first time qualitatively discussed the universal critical behavior of the the fourth order (kurtosis) of multiplicity fluctuations near the QCD critical point, which are realized by the coupling of particles with the order parameter σ field. The fluctuations of order parameter field σ(x) near a critical point can be described by the probability distributions as

$P[\sigma ]\sim \exp \{-\Omega [\sigma (x)]/T\}\mathrm{,}$ (49)

where Ω is the effective action functional for the field σ and can be expanded in the powers of σ,

$\Omega =\text{}{\displaystyle \int }\text{}{d}^{3}x\left[\frac{{(\nabla \sigma )}^2}{2}+\frac{m_{\sigma }^2}{2}{\sigma }^2+\frac{{\lambda }_3}{3}{\sigma }^3+\frac{{\lambda }_4}{4}{\sigma }^4+\dots \right]\mathrm{,}$ (50)

where mσ=1/ξ and the critical point is characterized by ξ → ∞. For the moments of the zero momentum mode ${\sigma }_V={\displaystyle \int }\sigma (x)d^{3}x$. Then, we have

$\langle {\sigma }_V^2\rangle =VT{\xi }^2\mathrm{,}$ (51) $\langle {\sigma }_V^3\rangle =2{\lambda }_{3}VT{\xi }^6\mathrm{,}$ (52) ${\langle {\sigma }_V^4\rangle }_c=6V{T}^3[2({\lambda }_3\xi {)}^2-{\lambda }_4]{\xi }^8\mathrm{,}$ (53)

where ${\langle {\sigma }_V^4\rangle }_c$ is the fourth order cumulants of the σ field. It is found that the higher order fluctuations are with higher power of the correlation length and diverge faster. If we introduce the coupling of the particles with the σ field, the fourth order cumulants of the particle multiplicity distributions can be obtained as

${\langle {(\delta N)}^4\rangle }_c=\langle N\rangle +{\langle {\sigma }_V^4\rangle }_c{\left(\frac{gd}{T}{\displaystyle {\int }_p}\frac{n_p}{{\gamma }_p}\right)}^4+\dots \mathrm{,}$ (54)

where np is the equilibrium distributions for a particle of a given mass, γp=(dEp/dm)^-1 is the relativistic gamma factor of a particle with momentum p and mass m, g is the coupling constant and d is the degeneracy factor. The mean value 〈N〉 in the r.h.s. of the Eq.(54) is the pure statistical contribution (Poisson).

Figure 8 left displays the sketch of QCD phase diagram with critical contributions to the σ field. When the chemical freeze-out lines (green dashed line) pass by the critical point from the crossover side, the probability distributions of the σ field change from gaussian to the double-peak non-gaussian distribution and the corresponding fourth order cumulant change from zero to negative (red region) and to positive (blue region). When this σ field couples with the particles, it leads to a non-monotonic energy dependence of the normalized fourth order cumulants of multiplicity distributions (ω₄=〈(δN)⁴〉c/〈N〉) along the chemical freeze-out line, as shown in the right of the Fig. 8, where the baseline is unity, the Poisson baseline. However, one has to keep in mind that here we only consider the critical point and statistical fluctuation contributions. Other dynamical effects in heavy-ion collisions, such as the effects of baryon number conservations, hadronic scattering and resonance decay, are not taken into account. Furthermore, the finite size and finite time effects, non-equilibrium memory effects are also important and need to be carefully studied.

Figure 8:Full-screen View

(Color online) (Top left) The sketch of the QCD phase diagram with sign of the fourth order cumullants of the σ field due to critical contributions [18]. The red region represents negative values and the blue region are positive values. The green dashed line is the chemical freeze-out lines in heavy-ion collisions. (Bottom left) The probability distributions of the σ field and the corresponding sign of the fourth order cumulants of the σ field distributions. (Right) The expected non-monotonic energy dependence for the normalized fourth order cumulant of multiplicity distributions (ω₄=〈(δ N)⁴〉c/〈N〉) when the chemical freeze-out line passes the critical region indicated in the left-top plot.

This critical point induced non-monotonic energy dependence of the fourth order cumulants along the chemical freeze-out line has been confirmed by many other model calculations, such as NJL [57-59], PQM [25, 26], chiral hydrodynamics [60, 61], and other calculations [62-65]. It indicates the σ field calculations capture the main feature of the critical point. However, it is still a crude model. Here, the σ field model only considers the critical fluctuations in static and infinite medium without taking account for the off-equilibrium effects in the dynamical expanding of the fireball created in heavy-ion collisions. Recently, a theoretical paper discussed critical fluctuations considering the off-equilibrium effects within Kibble-Zurek framework and observed a universal scaling of critical cumulants [66].

4.1.2 NJL Model

A QCD based effective model-the so called Nambu-Jona-Lasinio (NJL) model is also widely used to study the conserved charge fluctuations near the QCD critical point. In this model, quark and gluon are the basic degree of freedom. Although, there is no mechanism of the quark confinement implemented in the NJL model, it is still a simple and useful way to study the qualitative behavior of the susceptibility of the conserved charges near the QCD critical point. Here, we just show the results of two susceptibility ratios calculated from NJL model

$m_1(q)=\frac{{\chi }_3^q}{{\chi }_2^q}\mathrm{,}{m}_2(q)=\frac{{\chi }_4^q}{{\chi }_2^q}\mathrm{,}$ (55)

where q=B,Q,S, ${\chi }_n^q$ is the n^th order susceptibility. Figure 9 shows the sign of the m₁ and m₂ of baryon, charge and strangeness number. The red region is of positive values whereas the blue region represents negative values. The yellow regions represent the values of m₁ and m₂ are very close to zero. One may notice that the signals from baryon number fluctuations are stronger than from charge and strangeness number fluctuations. This is mainly due to the mass effects that the strange quark mass (ms) is much heavier than the mass of light quarks (mu,md). Figure 10 displays three colored hypothetical chemical freeze-out lines. The red solid freeze-out line is fitted to recent experimental data. The parametrized formula for obtaining the three curves are

Figure 9:Full-screen View

(Color online) The sign of the m₁ (top) and m₂ (bottom) of the baryon (B), charge (Q) and strangeness (S) number. The red region represents positive value while blue region represents negative value. The dashed line is the crossover line while the crosses denotes the first order phase transition boundary [57].

Figure 10:Full-screen View

(Color online) The red soild, blue dot-dashed and green dash lines represent three hypothetical freeze-out curves [57]. The black dashed line is the crossover line and the crosses denote the curve of the first order phase transition boundary. The triangles are experimental chemical freeze-out data.

$T({\mu }_B)=a-b{\mu }_B^2-c{\mu }_B^4\mathrm{,}$ (56)

where a=0.158 GeV, b=0.14 GeV^-1, and c=0.04 (solid), 0.08 (dot-dashed), 0.12 (dashed) GeV^-3. The relation between baryon chemical potential (μB) and collision energy can be parametrized as [67]

${\mu }_B(\sqrt{s})=\frac{1.477}{1+0.343\sqrt{s}}\mathrm{.}$ (57)

With the freeze-out curve and Eq. (57), we plot the m₁,m₂ of B,Q,S as a function of colliding energy in Fig. 11 along the three chemical freeze-out lines as shown in Fig. 10. The black dashed lines in Fig. 11 left are the results from the free quark gas model. When approaching the critical point at low energies, the NJL model predicts non-monotonic signal of the susceptibility ratios while for the free gas case all moments are close to 0. Furthermore, we can infer that m₂(B) should be a better probe of the critical behavior due to larger magnitude in signal and also the most important one, having sign changes from negative to positive with respect to collision energy than other cases. As we mentioned, since there has no quark confinement in NJL model, the baselines obtained from NJL model (away from critical point) are different from the ones from hadron resonance gas model, which is unity. One can also see that the behavior near QCD critical point is very much different from the results of weakly interacting quark gas. The behavior of these two quantities m₁(B) and m₂(B) at colliding energies at few GeV where experiments have not covered yet are of great importance as some other models predict opposite slope of these two quantities compared to the NJL prediction. Figure 11 right shows correlations between the m₂ and m₁ for baryon, charge, and strangeness, respectively. We can see that the m₂ and m₁ correlation along the three chemical freeze-out lines for baryon shows a closed loop with sign changes and looks like a banana shape. This is very different behavior comparing with the charge and strangeness sector.

Figure 11:Full-screen View

(Color online) (Left) The m₁ and m₂ of baryon (B), charge (Q) and strangeness (S) as a function of colliding energy along the three hypothetical freeze-out lines as plotted in Fig. 10. The black dashed lines are the results from a free quark gas model. (Right) The correlation plot m₂ versus m₁ for baryon (top), charge (middle) and strangeness (bottom) along three hypothetical freeze-out lines [57].

4.2 Baselines and Background Effects in Heavy-ion Collisions

In this section, we discuss the statistical baselines and some of the non-CP physics background effects for the fluctuations measurements in heavy-ion collisions. The discussion of thermal blurring, diffusion and resonance decay effects can be found in Ref. [68-70] and Ref. [71, 72, 63, 64], respectively.

4.2.2 Effects of Baryon Number Conservation and Nuclear Potential on Net-Proton (Baryon) Cumulants

The effects of baryon number conservations (BNS) and mean field potential are more and more important at low energies. To study those effects on the fluctuations of net-proton (baryon) number, the rapidity and transverse momentum dependence for the cumulants of the net-proton (baryon) multiplicity distributions in Au+Au collisions at $\sqrt{s_{\text{NN}}}=5\text{GeV}$ have been studied within a microscopic hadronic transport (JAM) model [75]. The simulations were done with two different modes, which are the mean field and the softening of equation of state (EoS) mode, respectively. The softening of EoS is simulated by introducing attractive orbits in the two-body scattering to realize a smaller pressure of the system. It was found that the mean field potential and softening of EoS have strong effects on the rapidity distributions (dN/dy) and the shape of the net-proton (baryon) distributions. By comparing the results from the two modes with the results from default cascade, one found that the net-proton (baryon) cumulants and ratios from the three modes have similar trends and show strong suppression with respect to unity, which is attributed to the effects of baryon number conservations [76, 51]. It means that the effects of mean field potential and softening of EoS might be not responsible for the observed strong enhancement in the most central (0-5%) Au+Au collisions at 7.7 GeV measured by the STAR experiment at RHIC.

Figures 12 shows cumulant ratios of net-proton (baryon) distributions in Au+Au collisions at $\sqrt{s_{\text{NN}}}=5\text{GeV}$ from JAM model. When increasing the rapidity acceptance (Δ y), the net-proton (baryon) cumulant ratios will decrease, reach a minimum and then increase, which is the typical effects of baryon number conservation [76, 51]. For different net-proton (baryon) cumulant ratios, the position of the minimum are different. It indicates the mean field potential and softening of EoS will not lead to large increase above unity for the net-proton (baryon) cumulants ratios. Instead, due to the baryon number conservation, large suppression for the fluctuations of net-proton (baryon) are observed. The rapidity dependence for the cumulants ratios calculated from the three modes are with the similar trend. It suggests that the observed similar trends obtained by JAM model without implementing critical physics are dominated by the effects of baryon number conservation. On the other hand, one observes that the net-baryon cumulant ratios show larger suppression with respect to unity than the net-proton and the higher order cumulant ratios also show larger suppression than the lower order. On the other hand, as the mean field potential implemented in the JAM model is momentum dependent, it is also important to study the momentum dependence for the cumulants of net-proton distributions. In Fig. 13, for different transverse momentum range, we plot the cumulant ratios of net-proton distributions as a function of rapidity window, which are calculated with the three different modes. The results computed from different modes are with the similar trends. When the p_T coverage is enlarged, the cumulant ratios are suppressed with respected to unity, the Poisson expectations. When the p_T range is small, the fluctuations are dominated by Poisson statistics and the cumulant ratios are very close to unity. Another study for the effect of mean field on baryon number fluctuations done with a Relativistic Mean Field (RMF) approach can be found in Ref. [77].

Figure 12:Full-screen View

(Color online) Rapidity dependence for the Cumulants ratios of net proton and net baryon multiplicity distributions in Au+Au collisions at $\sqrt{s_{\text{NN}}}=5\text{GeV}$ from JAM model computed in the three different modes [75]. In the left shows σ²/M (C₂/C₁) and C₃/C₁. The figure in the right shows Sσ (C₃/C₂) and κσ²(C₄/C₂). The dashed horizontal lines are with the value of unity.

Figure 13:Full-screen View

(Color online) Rapidity dependence for the cumulants ratios of net-proton distributions in Au+Au collisions at $\sqrt{s_{\text{NN}}}=5$ GeV from JAM model computed in the three different modes and various transverse momentum ranges [75]. From top to bottom are σ²/M (C₂/C₁), C₃/C₁, Sσ (C₃/C₂), and κσ²(C₄/C₂), respectively. The dashed horizontal lines are with the value of unity.

4.2.3 Net-Proton versus Net-Baryon Kurtosis from UrQMD and AMPT Model

The STAR experiment measures net-proton fluctuations instead of net-baryon fluctuations and one may want to know to what extend they can reflect the net-baryon fluctuations in heavy-ion collisions. Therefore, fig. 3 demonstrates the comparison between moments of net-proton and net-baryon distributions from AMPT [78] and UrQMD [79] model calculations. We can find that the $\kappa {\sigma }^2$ of net-baryon distributions are systematically lower than the net-proton results. The differences are even bigger for low energies than high energies. There are two possible effects for the difference between net-proton and net-baryon fluctuations, one is the non inclusion of neutrons in the net-proton fluctuations, and the other one is the nucleon isospin exchanging process due to Δ resonance formation via pπ and nπ interaction, the so called isospin randomization, which will modify the net-proton fluctuations after the chemical freezeout. A set of formulas have been derived to convert the measured net-proton cumulants to the net-baryon cumulants by taking into account the above two effects [80, 81]. The converting formulas for various order net-baryon cumulants can be written as

$C_1^{\text{net}-B}=2C_1^{\text{net}-p}\mathrm{,}$ (58) $C_2^{\text{net}-B}=4C_2^{\text{net}-p}-2C_1^{\text{tot}-p}\mathrm{,}$ (59) $C_3^{\text{net}-B}=8C_3^{\text{net}-p}-12(C_2^p-C_2^{\overline{p}})+6C_1^{\text{net}-p}\mathrm{,}$ (60) $C_4^{\text{net}-B}=16C_4^{\text{net}-p}+16C_3^{\text{tot}-p}-64(C_3^p+C_3^{\overline{p}})+48C_2^{\text{net}-p}+12C_2^{\text{tot}-p}-26C_1^{\text{tot}-p}\mathrm{,}$ (61)

where tot-p means proton number plus anti-proton number. The right side of Fig. 14 shows the net-baryon $\kappa {\sigma }^2$ (C₄/C₂) results, converted from the net-protons fluctuations. Within large uncertainties, they are consistent with the net-baryon results directly calculated with the AMPT model.

Figure 14:Full-screen View

(Color online) Energy dependence of $\kappa {\sigma }^2$ of net-proton and net-baryon distributions for 0-5% Au+Au collisions from the UrQMD (left) and the AMPT string melting model (right). The results marked as solid black stars are based on theoretical calculations using Asakawa and Kitazawa’s formula. The error calculation is based on the Bootstrap method.

4.2.4 Cumulants and Correlation Functions

Fluctuations and correlations are closely related to each other and they are two sides of coins. The various order cumulants can be expressed into the linear combinations of the multi-particle correlation functions [82, 83], which are directly related to the correlation length (ξ) of system. The multi-particle density are related to factorial moments as

$F_n=\langle N(N-1)\mathrm{...}(N-n+1)\rangle ={\displaystyle \int }d{p}_1\mathrm{...}d{p}_n\rho (p_1\mathrm{,...,}{p}_n)\mathrm{,}$ (62)

where Fn is the n^th order factorial moment and ρ (p₁,...,pn) is the n particle density distributions. The integral sums over the interested phase space. The generation function of the factorial cumulants is

$g(t)=\ln \langle {(1+t)}^N\rangle ={\displaystyle \sum _{k=1}^{\infty }}{c}_k\frac{t^k}{k!}\mathrm{,}$ (63) $c_k={\frac{{\partial }^{k}g(t)}{\partial t^k}|}_{t=0}\mathrm{,}$ (64)

where ck is the kth order factorial cumulant, N is the random variable. We have the relation between factorial moments and correlation function as:

$F_1={\displaystyle \int }dp\rho (p)=\langle N\rangle \mathrm{,}$ (65) $F_2={\displaystyle \int }d{p}_{1}d{p}_2\rho (p_1\mathrm{,}{p}_2)=F_1{}^2+c_2\mathrm{,}$ (66) $\begin{array}{l}{F}_3={\displaystyle \int }d{p}_{1}d{p}_{2}d{p}_3\rho (p_1\mathrm{,}{p}_2\mathrm{,}{p}_3)\\ =F_1{}^3+3c_2{F}_1+c_3\mathrm{,}\end{array}$ (67) $\begin{array}{l}{F}_4={\displaystyle \int }d{p}_{1}d{p}_{2}d{p}_{3}d{p}_4\rho (p_1\mathrm{,}{p}_2\mathrm{,}{p}_3\mathrm{,}{p}_4)\\ =F_1{}^4+6c_2{F}_1^2+4F_1{c}_3+3c_2^2+c_4\mathrm{.}\end{array}$ (68)

On the other hand, the relation between factorial cumulant (ck) and cumulants can be expressed as

$c_k={\displaystyle \sum _{i=0}^k}{s}_1(k\mathrm{,}i)C_i\mathrm{,}$ (69) $C_k={\displaystyle \sum _{i=0}^k}{s}_2(k\mathrm{,}i)c_i\mathrm{,}$ (70)

where the s₁ is the sterling number of the first kind, Ci is i^th order cumulant. Then, we have the following equations:

$c_1=C_1=\langle N\rangle \mathrm{,}$ (71) $c_2=C_2-\langle N\rangle \mathrm{,}$ (72) $c_3=C_3-3C_2+2\langle N\rangle \mathrm{,}$ (73) $c_4=C_4-6C_3+11C_2-6\langle N\rangle \mathrm{,}$ (74) $C_1=c_1=\langle N\rangle \mathrm{,}$ (75) $C_2=\langle N\rangle +c_2\mathrm{,}$ (76) $C_3=\langle N\rangle +3c_2+c_3\mathrm{,}$ (77) $C_4=\langle N\rangle +7c_2+6c_3+c_4\mathrm{.}$ (78)

It is well known that high order cumulants (Cn,n>2) are zero for the gaussian distribution and thus these are ideal probe of the Non-Gaussianity. For correlation function cn (n>1), they are zero for Poisson distributions, thus can be used to measure the deviation from Poisson fluctuations. If we define the correlation strength parameter ${\widehat{c}}_k$ as

${\widehat{c}}_k=\frac{c_k}{{\langle N\rangle }^k}\mathrm{,}$ (79)

where k=2,3,4...,n. The different order correlation strength parameter ${\widehat{c}}_k$ can reflect different physics process of the system. If the system consists of many independent sources, the correlation strength will be diluted and it scales with the multiplicity as

${\widehat{c}}_k\propto \frac{1}{{\langle N\rangle }^{k-1}}\mathrm{.}$ (80)

For example, it is the case that the A+A system is superposed by many p+p collisions. However, if the particle sources are strongly correlated with each other, which is the case near the critical point, then we have

${\widehat{c}}_k\propto const\mathrm{.}$ (81)

The long range correlation become dominated near the critical point and the cumulant are dominated by the highest order correlation function as: Ck≈ck ∝〈N〉k. If the thermal statistical fluctuations dominated in the system, we have Ck≈ck ∝〈N〉. Thus, to search for the critical point in heavy-ion collisions, it is also important to study the centrality, energy and the rapidity dependence of the multi-particle correlation functions. This is effective way to look for the pattern of long range correlations near the critical point in heavy-ion collisions and it is also very useful to study the contributions of the non-critical backgrounds, such as the baryon number conservations, resonance decay, and hadronic scattering.

Figure 15 shows the rapidity dependence of the proton cumulants and correlation functions in Au+Au collisions at $\sqrt{s_{{}_{\text{NN}}}}$ =5 GeV from JAM transport model calculation [84]. It can be found that the second and third order correlation functions show negative values when enlarging the rapidity acceptance. The large negative values of the second and third order proton correlation functions also leads to strong suppression of the fourth order proton cumulant. Those observations can be understood in terms of the baryon number conservations. As one can see in Fig. 15, it seems that the fourth order proton correlation function is consistent with zero, which is due to the absence of long range correlations in this model calculations. It means the baryon number conservations, which is a large background effect for searching for the critical point with fluctuations of conserved quantities in heavy-ion collisions, has negligible effects on the fourth order proton correlation functions. In other words, due to insensitive to the baryon number conservations, the fourth order proton correlation function is an ideal probe of the long range correlations induced by the critical point.

Figure 15:Full-screen View

(Color online) Proton cumulants (top panels) and correlation functions (bottom panels) as a function of mean proton number (〈Np〉) in Au+Au collisions at $\sqrt{s_{{}_{\text{NN}}}}$ =5 GeV from JAM model. The mean proton number is varied by changing the rapidity coverage.

5.1 Collision Geometry and Centrality Definition

Before introducing the background suppression methods, we would like to firstly discuss about the centrality definition used in the heavy-ion collisions. The definition of the collision centralities for two colliding nuclei is not unique and can be defined by different quantities. A commonly used quantity is the so called impact parameter b, defined as the distance between the geometrical centers of the colliding nuclei in the plane transverse to their direction. Other quantities, such as the number of participant nucleons, N_part and the number of binary collisions, N_coll, can be also used. Figure 16 shows a Glauber Monte Carlo event of Au+Au collision at $\sqrt{s_{{}_{\text{NN}}}}$ = 200 GeV with impact parameter b=6 fm [89]. The blue and red solid circles represent the participant nucleons from the two colliding gold nuclei. Figure 17 shows the average number of participant nucleons (〈N_part〉) and average number of binary collisions (〈N_coll〉) as a function of impact parameter b. One can see that there is no one-to-one correspond between N_part, N_coll and impact parameter b. Unfortunately, all of those geometrical variables can’t be directly measured in the heavy-ion collision experiment. Since the particle multiplicity can be easily measured and also can reflect the initial geometry of heavy-ion collision. The centrality in heavy-ion collisions is usually determined by a comparison between experimental measured particle multiplicity and Glauber Monte-Carlo simulations [89]. It is denoted as a percentage value (for e.g. 0-5%, 5-10%,...) for a collection of events to represent the fraction of the total cross section. Figure 18 illustrates how to define a collision centrality in heavy-ion collisions by comparing the particle multiplicities with Glauber Monte Carlo simulation and the correlation between the particle multiplicities and the Glauber calculated quantities b and N_part. However, the relation between measured particle multiplicities and collision geometry is not one-to-one correspondence and there are fluctuations in the particle multiplicity even for a fixed collision geometry.

Figure 16:Full-screen View

(Color online) Illustration of a Glauber Monte Carlo event for Au+Au at $\sqrt{s_{\text{NN}}}=200$ GeV with impact parameter b = 6 fm in the transverse plane (left panel) and along the beam axis (right panel) [89]. The nucleons are drawn with a radius $\sqrt{{\sigma }_{\text{inel}}^{\text{NN}}/\pi }/2$. Darker disks represent participating nucleons.

Figure 17:Full-screen View

(Color online) Average number of participants (〈N_part〉) and binary nucleon-nucleon collisions (〈N_coll〉) along with event-by-event fluctuation of these quantities in the Glauber Monte Carlo calculation as a function of the impact parameter b [89].

Figure 18:Full-screen View

(Color online) An illustrated example of the correlation of the total inclusive charged-particle multiplicity N_ch with Glauber-calculated quantities(b, N_part) [89]. The plotted distribution and various values are illustrative and not actual measurements..

5.2 Centrality Bin Width Correction

The centrality bin width effect is caused by the volume variation within a wide centrality bin and will cause an artificial centrality dependence for the fluctuation observables [85, 86]. The centrality bin width correction (CBWC) is to suppress the volume fluctuations effects in the event-by-event fluctuation analysis within finite centrality bin width. Experimentally, measurements are usually reported for a wide centrality bin (a range of particle multiplicity), such as 0-5%, 5-10%,...etc., to reduce statistical errors. We know that the smallest centrality bin is determined by a single value of particle multiplicity. To suppress the centrality bin width effect in a wide centrality bin, we calculate the cumulants (Cn) for each single particle multiplicity bin (N_ch). Then, the results reported for this wide centrality bin (N_ch) is to take the weighted average. The weight is the corresponding number of events in the particle multiplicity bin divided by the total events of the wide centrality bin. The method can be expressed as

where the nr is the number of events for multiplicity bin r and the corresponding weight for the multiplicity r, ${\omega }_r=n_r/{\displaystyle \sum _{r=N_1}^{N_2}}{n}_r$. N₁ and N₂ are the lowest and highest multiplicity values for one centrality bin. Once having the centrality bin width corrected cumulants via Eq.(82), we can calculate the various order cumulant ratios, for exampleκσ²=C₄/C₂ and Sσ=C₃/C₂, where the κ and S are kurtosis and skewness, respectively. The final statistical error of cumulants and cumulant ratios for wide centrality bin can be calculated by standard error propagation.

To demonstrate the centrality bin width effect and test the method of centrality bin width correction, we have calculated the cumulants of net-proton distributions in Au+Au collisions from UrQMD model in different ways. Figure 19 show the centrality dependence of the cumulant ratios (Sσ, κσ²) of net-proton multiplicity distributions for Au+Au collisions at $\sqrt{s_{{}_{\text{NN}}}}$ =7.7, 11.5, 19.6, 27, 39, 62.4, and 200 GeV from the UrQMD model calculations. The open circle and open cross in Fig. 19 represent the results obtained with and without applying the CBWC in the nine centralities (0-5%, 5-10%, 10-20%, 20-30%...70-80%), respectively. For the nine centralities, we clearly observe that the results with CBWC are very different from those without CBWC. This indicates the volume fluctuations in one wide centrality bin do have a significant impacts on the value of cumulants and the CBWC will make the values of the cumulants systematically lower by reducing the effects of volume fluctuations in one wide centrality bin. The solid circles show the centrality dependence for 32 finer centrality bins (0-2.5%, 2.5-5%, 5-7.5%...77.5%-80%) without CBWC. In the case of 32 centrality bin, due to the finer bin width, the centrality bin width effects are expected to be very small. Interestingly, we found that the results calculated from 32 centrality bins show good agreement with the results from nine centralities with CBWC. This further confirms the effectiveness of the CBWC method described above. On the other hand, we also tried to use the statistical errors (error) as weight to perform the CBWC by replacing the weight factor nr in Eq. (82) with 1/error² for each single multiplicity bin. The statistical error can be obtained by the Delta theorem and/or bootstrap methods at each multiplicity bin. It is found that the $S\sigma$ with CBWC using the statistical error as weight are consistent with the results with events number weighted CBWC, but not for κσ². It means the error weighted method can not be used for CBWC, which may be due to the statistical error is not only related to the number of events but also the cumulants itself.

Figure 19:Full-screen View

(Color online) The centrality dependence of the moments products Sσ (Left) and κσ² (Right)of net-proton multiplicity distributions for Au+Au collisions at $\sqrt{s_{{}_{\text{NN}}}}$ =7.7, 11.5, 19.6, 27, 39, 62.4, 200GeV in UrQMD model [86]. The solid dots represent the results calculated from 32 centrality bins.

5.3 Volume Fluctuations Effects

Volume fluctuations are long standing notorious background for the event-by-event fluctuation analysis in heavy-ion collisions [86, 90-95]. This is originated from that one cannot directly measure the collision centrality and/or initial collision geometry of the system of two nuclei. It is difficult to completely eliminated as it is usually convoluted with the real fluctuation signals. Consequently, this will lead to undesirable volume fluctuations in the event-by-event fluctuation analysis of particle multiplicity in heavy-ion collisions. The volume fluctuations will enhance the values of cumulants of the the event-by-event multiplicity distributions. However, the model calculations in the paper [95] conclude that the effects of volume fluctuations is too small to explain the large increase found in the preliminary result of 0-5% most central net-proton $\kappa {\sigma }^2$ in Au+Au collisions at $\sqrt{s_{{}_{\text{NN}}}}$ =7.7 GeV measured by the STAR experiment.

In the following, we will demonstrate the volume fluctuations in net-proton multiplicity fluctuations from Au+Au collisions by using UrQMD model simulations and discuss the method to suppress the volume fluctuations. To avoid auto-correlation, the centrality are defined with charged particle multiplicities by excluding the protons and anti-protons used in the analysis. The relation between measured particle multiplicity and impact parameter is not one-to-one correspond and there are fluctuations in the particle multiplicity even for a fixed impact parameter. Thus, we could obtain a finite resolution of initial collision geometry by using particle multiplicity to determine the centrality. As the N_part can reflect the initial geometry (volume) of the colliding nuclei, the σ²/M of N_part distributions can be regarded as the centrality resolution for a certain centrality definition. Figure20 shows the centrality dependence of σ²/M of number of participant nucleons (N_part) distributions for Au+Au collisions at $\sqrt{s_{{}_{\text{NN}}}}$ =7.7 and 200 GeV in UrQMD calculations with four different centrality definitions. The different centrality definitions are corresponding to the charged particles with four different η coverage (|η|<0.5,1.0,1.5,2.0). It shows that when we define the centrality with |η|<2, the σ²/M of N_part distributions for 7.7 and 200 GeV are similar. We can find that more particles are used in the centrality determination, the better centrality resolution and smaller fluctuation of the initial geometry (volume fluctuation) we obtain. Figure 21 shows the energy dependence of moment product (Sσ, κσ²) of net-proton multiplicity distributions in Au+Au collisions from UrQMD calculations for three different centralities (0-5%, 30-40%, 70-80%) with four different η ranges for centrality definitions. The centrality bin width corrections have been applied for all of the cases. Different centrality definitions could cause different results due to changing of magnitude of the volume fluctuations. By extending the η coverage of the charged particles used in centrality definition, we find the volume fluctuations are strongly suppressed. The κσ² (fourth order fluctuation) is more sensitive to the volume fluctuations than the $S\sigma$ (third order fluctuation). On the other hand, the volume fluctuations have much smaller effects in the most central collisions (0-5%) than in peripheral and mid-central collisions. With those studies, we conclude that having more particles in the centrality definition is an effective way to improve the centrality resolution and suppress the effects of volume fluctuations on the event-by-event fluctuations observables in heavy-ion collisions.

Figure 20:Full-screen View

(Color online) The centrality dependence of the σ² / M of N_part distributions for Au+Au collisions at $\sqrt{s_{{}_{\text{NN}}}}$ =7.7 and 200 GeV in UrQMD model [86]. Four different centrality definitions are corresponding to the charged particles with different η coverage (|η|<0.5,1.0,1.5,2.0).

Figure 21:Full-screen View

(Color online) The energy dependence of the moments products (Sσ, κσ²) of net-proton multiplicity distributions for Au+Au collisions at $\sqrt{s_{{}_{\text{NN}}}}$ =7.7, 11.5, 19.6, 27, 39, 62.4, 200 GeV in UrQMD model [86]. Four different centrality definitions are corresponding to the charged particles with different η coverage (|η|<0.5,1.0,1.5,2.0).

In principle, both the centrality bin width effects and centrality resolution effects are originated from volume fluctuations. The former is the volume variation within one wide centrality bin, and the latter is due to the initial volume fluctuations. These are two different effects and should be treated separately. The centrality bin width effects not only depend on the bin size but also dependent on the centrality resolutions (or the way to define the centrality). In this sense, these two effects are related with each other and both depend on the centrality definition. Figure 22 shows the $S\sigma$ and $\kappa {\sigma }^2$ of net-proton distributions in Au+Au collisions from UrQMD model calculations with centrality definitions at different η coverage (|η|< 0.5 and 2). The larger η coverage means more particles are included in the centrality definition and with better centrality resolution. Indeed, it shows that the results from wider η coverage centrality definition are get suppressed comparing with the results from narrower η coverage centrality definition. At fixed centrality definition, the results with CBWC are always smaller than the results without CBWC.

Figure 22:Full-screen View

(Color online) The centrality dependence of the moments products Sσ (Left) and κσ² (Right)of net-proton multiplicity distributions for Au+Au collisions at $\sqrt{s_{{}_{\text{NN}}}}$ =7.7, 11.5, 19.6, 27, 39, 62.4, 200 GeV in UrQMD model [86]. The open circles and squares are the results with CBWC and without CBWC at |η| < 0.5 for the centrality definition, respectively. The solid circles and squares are the results with at |η| < 2 in the centrality definition, respectively.

5.4 Auto-correlation Effects

The auto-correlation effect is a background effect in the fluctuation analysis and will suppress the magnitude of the signals. For example, in net-proton fluctuation analysis, to avoid the auto-correlation, we should exclude the corresponding protons and anti-protons from the centrality definition. For net-kaon fluctuations, we need to exclude K⁺ and K^- in the centrality definition. To illustrate this effects, we calculate the net-proton fluctuations in Au+Au collisions from UrQMD model with two different centrality definitions. One is using all charged particles and the other use the multiplicity of only charged kaon and pion to define the collision centrality. Figure 23 shows that for $S\sigma$ and $\kappa {\sigma }^2$ of net-proton distributions, the results with auto-correlation are smaller than the ones without auto-correlation. Meanwhile, the auto-correlation effects are stronger at lower energies. This is because the overlap fraction of proton/anti-protons used in the fluctuation analysis and in the centrality definition increase when decreasing the energies. In the data analysis, to avoid auto-correlation, we have to exclude the particles used in the fluctuation analysis from the centrality definition.

Figure 23:Full-screen View

(Color online) The centrality dependence of the moments products Sσ (Left) and κσ² (Right)of net-proton multiplicity distributions for Au+Au collisions at $\sqrt{s_{{}_{\text{NN}}}}$ =7.7, 11.5, 19.6, 27, 39, 62.4, 200 GeV in UrQMD model with the two different centrality definitions [86].

5.5 Efficiency Correction for Cumulants

The detector always have a finite particle detection efficiency. The observed event-by-event particle multiplicity distributions are the convolution between the original distributions and the efficiency response function. We need to correct this efficiency effect and a deconvolution operation is needed to recover the true fluctuations signals. However, it is not straightforward to get the efficiency corrected results for the cumulants of particle multiplicity distributions, especially for the higher order fluctuations.

It is well know that the detection efficiency response function is binomial distribution for a detector with good performance. Based on binomial efficiency response function, there has many discussions about the efficiency correction methods for moment analysis [87, 88]. Here, we provide a unified description of efficiency correction and error estimation for cumulants of multiplicity distributions [56]. The principle idea is to express the moments and cumulants in terms of the factorial moments, which can be easily corrected for efficiency effect. By knowing the covariance between factorial moments, we use the standard error propagation based on the Delta theorem in statistics to derive the error formulas for efficiency corrected cumulants. More important, this method can be also applied to the phase space dependent efficiency case, where the efficiency of proton or anti-proton are not constant within studied phase space. One needs to note that the efficiency correction and error estimation should be done for each single particle multiplicity bin in each centrality and just before the centrality bin width correction.

In the STAR experiment, the particle detection efficiency can be obtained from the so called Monte Carlo (MC) embedding techniques [96]. The Monte Carlo tracks are blended into real events at the raw data level. The tracks are propagated through the full simulation chain of the detector geometry with a realistic simulation of the detector response. The efficiency can be obtained by the ratio of matched MC tracks to input MC tracks. It contains the net effects of tracking efficiency, detector acceptance, decays, and interaction losses. For illustration purpose, we discuss the application of the efficiency correction on the net-proton fluctuation analysis in heavy-ion collisions. Experimentally, we measure net-proton number event-by-event wise, $n=n_p-n_{\overline{p}}$, which is proton number minus anti-proton number. The average value over the whole event ensemble is denoted by 〈n〉, where the single angle brackets are used to indicate ensemble average of an event-by-event distributions. For simplify, let us discuss constant efficiency case for (anti-)proton within the entire phase space. The probability distribution function of measured proton number np and anti-proton number $n_{\overline{p}}$ can be expressed as Ref. [87]:

where the $P(N_p\mathrm{,}{N}_{\overline{p}})$ is the original joint probability distribution of number of proton (Np) and anti-proton ($N_{\overline{p}}$), εp and $s{\epsilon }_{\overline{p}}$ are the efficiency of proton and anti-proton, respectively. To derive the efficiency correction formula for moments and cumulants, let us introduce the bivariate factorial moments:

With the Eqs. (83), (84) and (85), one can obtain a useful relation between the efficiency corrected and uncorrected factorial moments as

Then, the various order moments and cumulants can be expressed in terms of the factorial moments. Before deriving the formulas for the moments and cumulants of net-proton distributions, we need some mathematical relationships between moments, central moments, cumulants and factorial moments. Let us define a multivariate random vector 𝑿=(X₁,X₂,...,Xk)′ and a set of number 𝒓=(r₁,r₂,...,rk′. The multivariate moments, central moments and factorial moments can be written as

where E denotes the expectation value operator, and the L-Endmr(X), μr(X) and Fr(X), are multivariate moments, central moments, and factorial moments, respectively. Then, we have the relation between the moments and central moments by using bin omial expansions

where i=(i₁, i₂,..., ik)′. To get the relation between moments and factorial moments, one needs the Stirling numbers of the first (s₁(n,i)) and second kind (s₂(n,i)), which are defined as

where N, n and i are non-negative integer number. The recursion equations for the Stirling numbers of the first and second kind are

and