Search for the QCD Critical Point with Fluctuations of Conserved Quantities in Relativistic Heavy-Ion Collisions at RHIC : An Overview

TL;DR

This paper reviews experimental measurements of fluctuations in conserved quantities in heavy-ion collisions at RHIC to search for the QCD critical point, discussing data analysis, theoretical signatures, backgrounds, and future plans.

Contribution

It provides a comprehensive overview of fluctuation measurements, analysis techniques, and the implications for locating the QCD critical point in heavy-ion collisions.

Findings

Cumulant measurements show non-monotonic behavior near certain energies.

Analysis methods effectively suppress volume fluctuations and auto-correlations.

Future data from the second phase of RHIC's beam energy scan will enhance critical point search.

Abstract

Fluctuations of conserved quantities, such as baryon, electric charge and strangeness number, are sensitive observables in relativistic heavy-ion collisions to probe the QCD phase transition and search for the QCD critical point. In this paper, we review the experimental measurements of the cumulants (up to fourth order) of event-by-event net-proton (proxy for net-baryon), net-charge and net-kaon (proxy for net-strangeness) multiplicity distributions in Au+Au collisions at $\sqrt{s_{NN}}=7.7, 11.5, 14.5, 19.6, 27, 39, 62.4, 200$ GeV from the first phase of beam energy scan program at the Relativistic Heavy-Ion Collider (RHIC). We also summarize the data analysis methods of suppressing the volume fluctuations, auto-correlations and the unified description of efficiency correction and error estimation. Based on theoretical and model calculations, we will discuss the characteristic signatures of critical point as well as backgrounds for the fluctuation observables in heavy-ion collisions. The physics implications and the future second phase of the beam energy scan (2019-2020) at RHIC will be also discussed.

Tables

Table 1. Table 1: Locations of the QCD critical point from Lattice QCD and DSE, respectively.

	Lattice			DSE
(${\mu }_B^E,T^E$) MeV	I qcp_zoltan	II qcp_Gupta2 ; Gupta_Lattice2016	III Karsch_QM2015 ; Karsch_CPOD2016 ; Karsch_SQM2016 ; Karsch_INT_2016	I xin2014quark	II shi2014locate	III fischer2014phase
	(360,162)	(285,155)	${\mu }_B^E/T^E$>2	(372,129)	(405,127)	(504,115)

Equations

$$\frac{P}{T^4}=\frac{1}{V{T}^3}\ln [Z(V,T,{\mu }_B,{\mu }_Q,{\mu }_S)]$$

$${\chi }_{ijk}^{BQS}=\frac{{\partial }^{(i+j+k)}[P/T^4]}{\partial {\hat{\mu }}_B^i\partial {\hat{\mu }}_Q^j\partial {\hat{\mu }}_S^k}$$

$$\begin{array}{rl}{C}_{ijk}^{BQS}& =\frac{{\partial }^{(i+j+k)}\ln [Z(V,T,{\mu }_B,{\mu }_Q,{\mu }_S)]}{\partial {\hat{\mu }}_B^i\partial {\hat{\mu }}_Q^j\partial {\hat{\mu }}_S^k}\\ & =V{T}^3{\chi }_{ijk}^{BQS}(T,{\mu }_B,{\mu }_Q,{\mu }_S)\end{array}$$

$$M_q$$

$$C_3^q$$

$$C_5^q$$

$$C_6^q$$

$$S_q$$

$${\kappa }_q$$

$$\frac{{\sigma }_q^2}{M_q}$$

$${\kappa }_q{\sigma }_q^2$$

$$\frac{P}{T^4}$$

$${\chi }_{2n}^B$$

$${\chi }_{2n-1}^B$$

$$\frac{C_{\mathrm{even}}^B}{C_{\mathrm{even}}^B}$$

$$\frac{C_{\mathrm{odd}}^B}{C_{\mathrm{even}}^B}$$

$$\frac{M_B}{{\sigma }_B^2}$$

$${\kappa }_B{\sigma }_B^2$$

$$\begin{array}{rl}& \frac{P(T,{\mu }_B)-P(T,0)}{T^4}=\frac{1}{2}{\chi }_2^B(T){\left(\frac{{\mu }_B}{T}\right)}^2\times \\ & \left[1+\frac{1}{12}\frac{{\chi }_4^B(T)}{{\chi }_2^B(T)}{\left(\frac{{\mu }_B}{T}\right)}^2+\frac{1}{360}\frac{{\chi }_6^B(T)}{{\chi }_2^B(T)}{\left(\frac{{\mu }_B}{T}\right)}^4\right]+\mathcal{O}({\mu }_B^8)\end{array}$$

$${\chi }_n^B(T,{\mu }_B)$$

$${\chi }_n^Q(T,{\mu }_B)$$

$${\chi }_n^S(T,{\mu }_B)$$

$${\chi }_{2n-1}^B(T,{\mu }_B)$$

$${\chi }_{2n}^B(T,{\mu }_B)$$

$$L_n=\frac{1}{6}\frac{{\chi }_{2n+2}^B(T)}{{\chi }_{2n}^B(T)}{\left(\frac{{\mu }_B}{T}\right)}^2$$

$${\chi }_{2n-1}^B(T,{\mu }_B)$$

$${\chi }_{2n}^B(T,{\mu }_B)$$

$$\frac{{\chi }_{2n}^B(T,{\mu }_B)}{{\chi }_{2n-1}^B(T,{\mu }_B)}$$

$$\frac{{\chi }_{2n+1}^B(T,{\mu }_B)}{{\chi }_{2n}^B(T,{\mu }_B)}$$

$$\frac{{\chi }_{2n+2}^B(T,{\mu }_B)}{{\chi }_{2n}^B(T,{\mu }_B)}$$

$$\frac{{\chi }_{2n+1}^B(T,{\mu }_B)}{{\chi }_{2n-1}^B(T,{\mu }_B)}$$

$$\frac{{\chi }_{2n}^B(T,{\mu }_B)}{{\chi }_{2n-1}^B(T,{\mu }_B)}$$

$$\frac{{\chi }_{2n+1}^B(T,{\mu }_B)}{{\chi }_{2n}^B(T,{\mu }_B)}$$

$$\frac{{\chi }_{2n+2}^B(T,{\mu }_B)}{{\chi }_{2n}^B(T,{\mu }_B)}$$

$$\frac{{\chi }_{2n+1}^B(T,{\mu }_B)}{{\chi }_{2n-1}^B(T,{\mu }_B)}$$

$$\begin{array}{rl}& \frac{{\chi }_{2n+2}^B(T,{\mu }_B)}{{\chi }_{2n}^B(T,{\mu }_B)}-\frac{{\chi }_{2n+1}^B(T,{\mu }_B)}{{\chi }_{2n-1}^B(T,{\mu }_B)}=\\ & \frac{1}{3}\left[\frac{{\chi }_{2n+4}^B(T)}{{\chi }_{2n}^B(T)}-{\left(\frac{{\chi }_{2n+2}^B(T)}{{\chi }_{2n}^B(T)}\right)}^2\right]{\left(\frac{M_B}{{\sigma }_2^B}\right)}^2\end{array}$$

$$r_n(T)=\frac{1}{3}\left[\frac{{\chi }_{2n+4}^B(T)}{{\chi }_{2n}^B(T)}-{\left(\frac{{\chi }_{2n+2}^B(T)}{{\chi }_{2n}^B(T)}\right)}^2\right]$$

$$r_n(T)=\left(\frac{{\chi }_{2n+2}^B(T,{\mu }_B)}{{\chi }_{2n}^B(T,{\mu }_B)}-\frac{{\chi }_{2n+1}^B(T,{\mu }_B)}{{\chi }_{2n-1}^B(T,{\mu }_B)}\right){\left(\frac{{\sigma }_B^2}{M_B}\right)}^2$$

$$\frac{{\sigma }_B^2}{M_B}(T,{\mu }_B)$$

Full text

QCD critical point, Fluctuations and correlations, Relativistic heavy-ion collisions, Conserved charges

• I)

Introduction

• II)

The QCD Critical Point

• III)

Fluctuations and Correlations

• III.1)

Hadron Resonance Gas Model

• III.2)

Lattice QCD

• IV)

Experimental Observables

• IV.1)

Fluctuations Signature near the QCD Critical Point

• IV.1.1)

$\sigma$ Field Model

• IV.1.2)

NJL Model

• IV.2)

Baselines and Background Effects for Net-Particle Cumulants

• IV.2.1)

Expectations from Poisson, Binomial and Negative Binomial Statistics

• IV.2.2)

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

• IV.2.3)

Fluctuations of Net-Proton (Baryon) from UrQMD and AMPT

• IV.2.4)

Cumulants and Correlation Functions

• V)

Data Analysis Methods

• V.1)

Collision Geometry and Centrality Definition

• V.2)

Centrality Bin Width Correction

• V.3)

Volume Fluctuations Effects

• V.4)

Auto Correlation Effects

• V.5)

Efficiency Corrections for Cumulants

• V.6)

Error Estimations for Efficiency Corrected Cumulants

• VI)

Experimental Results

• VII)

Beam Energy Scan Phase-II and STAR Detector Upgrades

• VIII)

Summary

I Introduction

A major uncertainty in our understanding of strongly interacting nuclear matter is the so called Quantum Chromodynamics (QCD) phase structure and the possible existence of a critical point in the QCD phase diagram, located at high temperature and non-zero baryon chemical potential. It is one of the main goals of the Beam Energy Scan (BES) program at the Relativistic Heavy-Ion Collider (RHIC) bes ; BESII_WhitePaper , which is located at the Brookhaven National Laboratory (BNL), US. This also serves as a main motivation for the research programs at the future accelerator facilities FAIR in Darmstadt and NICA in Dubna. As shown in Fig.1, the conjectured QCD phase diagram, it can be displayed in the two dimensional phase diagram (temperature, $T$ vs. baryon chemical potential, $\mu_{B}$). Finite temperature Lattice QCD calculations has shown that at zero $\mu_{B}$ ($\mu_{B}=0$) region, it is a crossover transition between hadronic phase and quark-gluon plasma (QGP) phase crossover . At large $\mu_{B}$ region, the QCD based models predicted that the phase transition is of the first order de2002qcd ; endrHodi2011qcd and there should exist a so called QCD Critical Point (CP) as the endpoint of the first order phase boundary QCP_Prediction ; fodor2004critical . Due to sign problem at finite $\mu_{B}$ region, it is difficult to precisely determine the location of the CP or even its existence gavai2015qcd . Experimental confirmation of the existence of the CP will be an excellent verification of QCD theory in the non-perturbative region and a milestone of exploring the QCD phase structure. Please note that the first-order phase boundary, the critical point and the smooth crossover are closely related thermodynamically. For example, if the smooth crossover and the first-order phase boundary exist, there must be a critical point at the end of the first-order line. To some extent, the burden is on the experimental side who should determine the location of the QCD critical point or the first-order phase boundary. To access a broad region of the QCD phase diagram, experimentalists vary the temperature ($T$) and baryon chemical potential ($\mu_{B}$) of the nuclear matter created in heavy-ion collisions bes by tuning the colliding energies of two nuclei. It is expected that fluctuations of conserved charges yield information on the phase structure of QCD matter stephanov1998signatures ; stephanov1999event ; jeon1999fluctuations ; asakawa2000fluctuation ; koch2005baryon ; ejiri2006hadronic ; ratioCumulant , provided the freeze-out is sufficiently close to the phase boundary. These conserved quantities have been long time predicted to be sensitive to the correlation length qcp_signal ; ratioCumulant ; Neg_Kurtosis ; Hatta and directly connected to the susceptibilities computed in the first principle Lattice QCD calculations science ; ding2015thermodynamics ; Lattice ; MCheng2009 ; bazavov2012fluctuations ; bazavov2012freeze ; BFriman_EPJC ; 2014_Bengt_flu ; chiral_HRG ; kenji . Consequently, the analysis of event-by-event fluctuations of the net baryon number ($B$), electric charge ($Q$), and strangeness ($S$), in particular their higher order cumulants, play a central role in the efforts to reveal the thermodynamics of the matter created in heavy-ion collisions at RHIC and LHC. Thus, it can serve as a powerful observables to study the phase transition and search for the CP in heavy-ion collisions Misha2009_PRL ; Asakawa ,.

During the first phase of the RHIC BES (2010 to 2014), the STAR experiment has measured the cumulants (up to the fourth order) of net-proton (proton minus anti-proton number, proxy of net-baryon Hatta ) , net-charge and net-kaon multiplicity distributions in Au+Au collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}$= 7.7, 11.5, 14.5, 19.6, 27, 39, 62.4 and 200 GeV. In those energies, the data of 14.5 GeV is taken in the year 2014, 19.6, 27 GeV are taken in the year 2011, and the other energies are collected in the years 2010. In this paper, we will review the experimental results on fluctuations of conserved quantities from the BES data measured by the STAR and PHENIX experiments. The corresponding physics implications will be also discussed.

II The QCD Critical Point

A critical point is the end point of the first order phase transition boundary in the phase diagram, at which, the phase transition is of the second order and one cannot distinguish difference between the two phases. For e.g., in the liquid-gas phase transition of water, one cannot distinguish vapor and liquid of the water when the temperature is above the critical temperature $T_{c}$ (373.946 ℃). In equilibrated matter in the vicinity of a critical point, various thermodynamic quantities exhibit large critical fluctuations, which in laboratory systems give rise to e.g. critical opalescence. The critical phenomena (critical opalescence) is discovered by Baron Charles Cagniard de la Tour in 1822 in the study of the liquid-gas phase transition for the mixtures of alcohol and water berche2009critical . The term ”critical point” is firstly named by Thomas Andrews in 1869 andrews1869bakerian when he studied the liquid-gas transition in carbon dioxide ($\mathrm{CO_{2}}$), the critical temperature is about 31 ℃. When the thermodynamic condition of system is approaching the critical point, the correlation length of system will diverge. The divergency of the correlation length ($\xi$) is one of the most important characteristic feature of the critical point and it is also related to the divergency of the specific heat ($\mathrm{C_{v}}$), susceptibility ($\chi$), compressibility ($\kappa$) and critical opalescence. In Fig. 2, it demonstrates the well-known critical opalescence, the visible cloudy phenomena near the critical point of liquid-gas phase transition. When the lights are passing through the $\mathrm{CO_{2}}$ near the critical point, the light will undergo large scattering due to its wavelength is comparable to the length scale (correlation length) of the density fluctuations in the phase transition of the liquid-gas system.

Those critical behaviors can be described by power law divergence with a set of critical exponents. The critical exponents of the critical point for various systems with same symmetry and dimension belong to the same universality class. Due to self-similarity and scaling properties of the critical point, those critical exponents can be precisely calculated by the renormalization group theory wilson1974renormalization . Another important feature of the critical point is the so called finite size effect, which is originated from that the correlation length is comparable with the size of system and the system size limits the growing of the correlation length. This leads to an observable effects when one varies the system size.

The phase diagram of water is shown in Fig. 3 braun2009colloquium . It can be found that the phase structure of water are very rich, which is the emergent properties of quantum electrodynamics (QED). Due to the easily realized phase transition conditions, the water phase diagram are precisely known. On the other hand, the phase structure of the hot and dense nuclear matter, which is governed by the strong force described by the QCD theory, is rarely known to us. Thus, it is very important to explore the QCD phase structure and search for the QCD critical point theoretically and experimentally. From theoretical side, it is still very difficult to precisely determine the location of the critical point due to its non-perturbative feature. The QCD based models, such as NJL, PNJL, PQM, have given many results of the location of the QCD critical point, which are summarized in the reference location . The locations of the QCD critical point obtained from the first principle Lattice QCD and Dyson-schwinger equation (DSE) calculations are summarized in the table 1. One can see that the baryon chemical potential ($\mu_{B}^{E}$) of the QCD critical point are ranging from 266 to 504 MeV, the critical temperature is from 115 to 162 MeV. There still has big difference between the results from different methods and groups. Experimentally, we aim to search for the critical point with the strongly interacting QCD matter created in the relativistic heavy-ion collisions. It is very challenging due to the following reasons: 1. The hot dense medium created in the heavy-ion collisions are not static but expanding rapidly. Thus, the correlation length of the system is not only limited by the size of the system, but also by the finite expansion time and it is predicted to be 2-3 $fm$ by assuming the existence of a critical point correlationlength . One has to consider finite time and finite size effects in order to determine the exact location of the critical point. 2. What’s the sensitive observables and what’s the smoking gun signature of the QCD critical point in heavy-ion collisions. 3. One has to understand the non-critical contributions to the experimental observables and the signal to background ratio should not be too small. 4. One needs that the freeze-out thermodynamic conditions of the QCD matter created in heavy-ion collisions should be close enough to the phase boundary that the phase transition signals weren’t washed out after the expansion.

Due to the difficulties and challenges discussed above, we should set up good strategies to search for the QCD critical point. Firstly, we need to have good quality experimental data of heavy-ion collisions at a wide range of energies. This allows us to scan a broad region of the QCD phase diagram. Then, we use sensitive observables to find the smoking gun signatures and confirm the existence of the QCD critical point before determining its location. In order to extract critical signature and understand the background contributions, careful modelling of the critical phenomena and dynamical evolution of the heavy-ion collisions are needed. It requires close collaboration between theorists and experimentalists. If the QCD critical point is given and hidden in nature, we will finally discover it and put a permanent landmark in the phase diagram of the strongly interacting nuclear matter.

III Fluctuations and Correlations

Fluctuations and correlations have long been considered to be sensitive observables in heavy-ion collisions to explore the phase structure of the strongly interacting QCD matter jeon2004event ; koch2010hadronic ; asakawa2000fluctuation . They have a well defined physical interpretation for a system in thermal equilibrium and can provide essential information about the effective degrees of freedom. The well known phenomenon of critical opalescence is a result of fluctuations at all length scales due to a second order phase transition. The most efficient way to study the fluctuations of a system created in a heavy-ion collision is to measure an observable on the event-by-event basis and the fluctuations are studied over the ensemble of the events. In strong interaction, the net number of charges in a closed system is conserved. The magnitude of these fluctuations in a grand canonical ensemble at finite temperature are distinctly different in the hadronic and quark gluon plasma phases. Event-by-event fluctuation and correlation of the conserved charges is one of the observables to study the properties of the QCD medium created in relativistic heavy-ion collisions. Although these observables are hadronic ones, it is believed that they can reflect the thermal property in the early stage. A system in thermal equilibrium (for a grand-canonical ensemble) can be characterized by its dimensionless pressure, which is the logarithm of the QCD partition function ding2015thermodynamics :

[TABLE]

where $V$ and $T$ are the system volume and temperature. The $\mu_{B}$, $\mu_{Q}$ and ${\mu_{S}}$ are baryon, charge and strangeness chemical potential, respectively. The equation of state is very different for thermodynamical system with different degree of freedom and interactions. The susceptibility of the conserved charges ($B,Q,S$) are defined as the derivative of the dimensionless pressure with respected to the reduced chemical potential.

[TABLE]

where $\hat{\mu_{q}}=\mu_{q}/T,q=B,Q,S$. The cumulants of these conserved quantities ($B,Q,S$) distributions are connected to the corresponding susceptibilities by

[TABLE]

where the $C_{ijk}^{BQS}$ denotes both diagonal and off-diagonal cumulants of conserved quantities ($B,Q,S$) ($i,j,k=1,2,3,4...n$). Experimentally, we construct the ratios of cumulants as the experimental observables, which cancel the volume dependent and can be directly compared with the ratios of susceptibilities from theoretical calculations. To obtain the ratio of cumulants, we firstly introduce various order cumulants up to sixth order and their relations to the central moments as:

[TABLE]

where $M_{q}$, $\sigma_{2}^{q}$ are the mean and variance, respectively. The $C_{n}^{q}(n=2,3,4,...)$ are the $n^{th}$ order cumulants with $q=B,Q,S$ and $\delta N_{q}=N_{q}-<N_{q}>$. We didn’t consider the correlations between different conserved charges.

On the other hand, we introduce two well known statistic quantities, the so called $skewness(S)$ and $kurtosis(\kappa)$. In statistics, those two quantities can be used to describe the shape of distributions and they are defined as :

[TABLE]

For gaussian distribution, both of the two quantities are equal to zero. Thus, they are widely used to quantify the non-gaussianity. With above definition of the mean, variance, skewness, kurtosis and various order cumulants, we can have the following relations:

[TABLE]

Those equations connect the experimental measurements ($l.h.s.$) and theoretically calculations ($r.h.s.$). In the following, we will discuss the results calculated from the Hadron Resonance Gas Model and Lattice QCD.

III.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 HRG_Karsch ; HRG_fjh ; HRG_FJH2 :

[TABLE]

where $g_{i}$ is the degeneracy factor for hadrons of mass $m_{i}$, and $\hat{\mu}_{q}\equiv\frac{\mu_{q}}{T}$, with $q=B$, $S$, $Q$ denote the net-baryon number, net-strangeness and the net-charge, and $\mu_{B}$, $\mu_{S}$, $\mu_{Q}$ are the corresponding chemical potentials respectively. The $K_{2}(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:

[TABLE]

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

[TABLE]

Based on the Eq.(19) and (20), we obtain:

[TABLE]

where $n=1,2,3...$, $\mu_{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. Fig.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 $\chi_{2}/\chi_{1}$ or $\chi_{3}/\chi_{2}$ of charge and strangeness show strong energy dependence whereas the $\chi_{4}/\chi_{2}$ has small variation with energies. Due to the contributions from the multi-charge states $Q=2$ or $S=2,3$, the charge and strangeness $\chi_{4}/\chi_{2}$ deviate from unity.

III.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 wilson1974confinement . 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.

Fig. 5 shows the results of QCD equation of state (the trace anomaly, the pressure and the entropy density) from two independent groups: HotQCD 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^{4}$) at finite $\mu_{B}$ region can be calculated by using the Taylor expansion techniques. By putting the $\mu_{Q}$=$\mu_{S}$=0, we can expand the pressure ($P/T^{4}$) into finite $\mu_{B}$ as C6_LatticeQCD ; ding2015thermodynamics ; Karsch_INT_2016 :

[TABLE]

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 $\mu_{B}=\mu_{Q}=\mu_{S}=0$. Fig.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. Fig. 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 $\mu_{B}/T<2$. It means the Taylor expansion up to NNLO order is under control with $\mu_{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 ($\mu_{B}$) at $\mu_{Q}=\mu_{S}=0$ C6_LatticeQCD ; ding2015thermodynamics ; Karsch_INT_2016 :

[TABLE]

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:

[TABLE]

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

[TABLE]

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

[TABLE]

Then, we can express the baryon number susceptibilities ratios as:

[TABLE]

If we consider $L_{n}<<1$, the $r.h.s.$ of the Eq.(32) to (35) can be simplified as :

[TABLE]

Based on the Eq.(36), (38) and (39), we have:

[TABLE]

where we use a leading order approximation $\frac{\sigma^{B}_{2}}{M_{B}}(T,\mu_{B})=\frac{T}{\mu_{B}}$. We define a temperature dependent quantity $r_{n}(T)$:

[TABLE]

Then, deriving from Eq. (40) and (41), we get:

[TABLE]

For the lowest order with $n=1$, we obtain:

[TABLE]

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:

[TABLE]

where

[TABLE]

In above Taylor expansions of the baryon number susceptibilities in Lattice QCD, we always assume $\mu_{Q}=\mu_{S}=0$ and expand up to next to leading order. For more realistic, one need to consider the case $\mu_{Q}\neq\mu_{S}\neq 0$. In order to compare with experimental data, we need additional constrains. For eg., strangeness neutrality ($N_{S}$=0) and baryon to charge number ratios equals to 0.4 ($N_{Q}/N_{B}$=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 Karsch_SQM2016 ; Karsch_INT_2016 .

IV Experimental Observables

Event-by-event particle multiplicity fluctuations can be characterized by the cumulants of the event-by-event multiplicity distributions. It can be calculated as

[TABLE]

where $N$ is particle or net-particle number measured on the event-by-event bias and the $\left<N\right>$ is average over entire event ensemble, $\delta N=N-\left<N\right>$. With the definition of cumulants, we can also define mean ($M$), variance ($\sigma^{2}$), skewness ($S$) and kurtosis ($\kappa$) as:

[TABLE]

In addition, the moments product $\kappa\sigma^{2}$ and $S\sigma$ can be expressed in terms of the ratios of cumulants:

[TABLE]

The ratios of cumulants are independent on system volume. The statistical errors of those cumulants and cumulants ratios are estimated by the Delta theorem Delta_theory ; Unified_Errors . In general, the statistical errors strongly depend on the shape of the distributions, especially the width. For gaussian distributions, the statistical errors of cumulants ($C_{n}$) can be approximated as $error(C_{n})\propto\sigma^{n}/(\sqrt{N}\epsilon^{n})$, where $\sigma$ is the measured width of the distribution, $N$ represents the number of events and $\epsilon$ is the particle detection efficiency. Theoretically, conserved charge fluctuations are sensitive to the correlation length ($\xi$) of system, which is about 2-3 $fm$ near the QCD critical point in heavy-ion collisions. The fourth order cumulant are proportional to seventh power of the correlation length as $C_{4}\propto\xi^{7}$. Experimentally, the various order cumulants of net-proton (proton number minus anti-proton number), net-charge and net-kaon multiplicity distributions are measured with the data of the beam energy scan at RHIC to search for the signature of the QCD critical point. Since the STAR detector at RHIC can not measure the neutron at mid-rapidity, the net-proton fluctuations are used to approximate the fluctuations of net-baryons. And the net-kaon multiplicity fluctuations also is used to approximate to the fluctuations of the net-strangeness. To search for the CP in heavy-ion collisions, the event-by-event multiplicity fluctuations are not necessary conserved quantities, for eg., the proton number fluctuations itself can reflect the singularity of the critical point and can be directly used to search for the CP. However, there are two advantages in using the conserved quantities, one is that it can be directly connected to the susceptibilities of system, which can be computed in first principle Lattice QCD, the other one is that due to the dynamical expansion of the QCD medium created in heavy-ion collisions, the signal of conserved quantities will not be easily washed out by the diffusion process and can be preserved in the final state Asakawa .

In the following, we will discuss the fluctuation signatures of the QCD critical point from various theoretical calculations, such as $\sigma$ field and NJL model. Finally, we will also show the effects of nuclear potential and baryon number conservations on the cumulants of net-proton (baryon) distributions.

IV.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.

IV.1.1 $\sigma$ Field Model

One of the most important calculations is done with the $\sigma$ field model Neg_Kurtosis . 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 $\sigma$ field. The fluctuations of order parameter field $\sigma(x)$ near a critical point can be described by the probability distributions as:

[TABLE]

where $\Omega$ is the effective action functional for the field $\sigma$ and can be expanded in the powers of $\sigma$:

[TABLE]

where $m_{\sigma}=1/\xi$ and the critical point is characterized by $\xi\to\infty$. For the moments of the zero momentum mode ${\sigma_{V}}=\int{\sigma(x){d^{3}}}x$. Then, we have

[TABLE]

where $\langle\sigma_{V}^{4}\rangle_{c}$ is the fourth order cumulants of the $\sigma$ 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 $\sigma$ field, the fourth order cumulants of the particle multiplicity distributions can be obtained as:

[TABLE]

where $n_{p}$ is the equilibrium distributions for a particle of a given mass, $\gamma_{p}=(dE_{p}/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 $\langle N\rangle$ in the $r.h.s.$ of the Eq.(58) is the pure statistical contribution (Poisson).

Figure 8 left displays the sketch of QCD phase diagram with critical contributions to the $\sigma$ field. When the chemical freeze-out lines (green dashed line) pass by the critical point from the crossover side, the probability distributions of the $\sigma$ 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 $\sigma$ field couples with the particles, it leads to a non-monotonic energy dependence of the normalized fourth order cumulants of multiplicity distributions ($\omega_{4}=\langle(\delta N)^{4}\rangle_{c}/\langle N\rangle$) 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.

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 2015_JianDeng_fluctuation ; Wenkai_NJL , PQM BFriman_EPJC ; 2014_Bengt_flu , chiral hydrodynamics Nahrgang_EPJA_2016 ; herold2016dynamical and other calculations VDW_flu ; Marcus_2016 ; Mukherjee_2016 . It indicates the $\sigma$ field calculations capture the main feature of the critical point. However, it is still a crude model. Here, the $\sigma$ 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 Swagato_2016_PRL .

IV.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 :

[TABLE]

where $q=B,Q,S$, $\chi^{q}_{n}$ is the $n^{th}$ order susceptibility. Fig. 9 shows the sign of the $m_{1}$ and $m_{2}$ of baryon, charge and strangeness number. The red region are of positive values wheres the blue region represents negative values. The yellow regions represent the values of $m_{1}$ and $m_{2}$ 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 ($m_{s}$) is much heavier than the mass of light quarks ($m_{u},m_{d}$). Fig.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

[TABLE]

where $a=0.158$ GeV, $b=0.14$ $\mathrm{GeV^{-1}}$, and $c=0.04$ (solid), $0.08$ (dot-dashed), $0.12$ (dashed) $\mathrm{GeV^{-3}}$. The relation between baryon chemical potential ($\mu_{B}$) and collision energy can be parametrized as Begun_SQM2016 :

[TABLE]

With the freeze-out curve and Eq. (61), we plot the $m_{1},m_{2}$ 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 [math]. Furthermore, we can infer that $m_{2}(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_{1}(B)$ and $m_{2}(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. Fig. 11 right shows correlations between the $m_{2}$ and $m_{1}$ for baryon, charge and strangeness, respectively. We can see that the $m_{2}$ and $m_{1}$ 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.

IV.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 Asakawa_review ; Asakawa_thermal ; 2015_Kitazawa_Rapidity and HRG_baseline ; HRG_Nahrgang ; Marcus_2016 , respectively.

IV.2.1 Expectations from Poisson, Binomial and Negative Binomial Statistics

In the following, we discuss some expectations for cumulants of net-proton multiplicity distributions from some basic distributions QM2014_baseline .

Poisson Distributions : If the particle and anti-particle are independently distributed as Poissonian distributions. Then the net-proton multiplicity will follow the Skellam distribution, which is expressed as:

$P(N)={(\frac{{{M_{p}}}}{{{M_{\overline{p}}}}})^{N/2}}{I_{N}}(2\sqrt{{M_{p}}{M_{\overline{p}}}})\exp[-({M_{p}}+{M_{\overline{p}}})],$ where the $N$ is the net-proton number, $I_{N}(x)$ is a modified Bessel function, $M_{p}$ and $M_{\overline{p}}$ are the mean number of particles and anti-particles, as shown in Fig. 1. The various order cumulants ($C_{n}$) are closely connected with the moments, e.g., $C_{1}=\langle N\rangle=M,C_{2}=\langle(\delta N)^{2}\rangle=\sigma^{2},C_{3}=\langle(\delta N)^{3}\rangle=S\sigma^{3},C_{4}=\langle(\delta N)^{4}\rangle-3\langle(\delta N)^{2}\rangle^{2}=\kappa\sigma^{4}$, where the $\delta N$ =$N-\langle N\rangle$, the $\sigma^{2}$, $S$ and $\kappa$ are variance, skewness and kurtosis, respectively. Then we construct, $S\sigma={C_{3}}/{C_{2}}=({M_{p}}-{M_{\overline{p}}})/({M_{p}}+{M_{\overline{p}}})$ and $\kappa{\sigma^{2}}={C_{4}}/{C_{2}}=1$, which provides the Poisson expectations for the various order cumulants/moments of net-particle distributions. The only input parameters of the Poisson baseline for cumulants of net-particle distributions are the mean values of the particle and anti-particle distributions. 2. 2.

Binomial and Negative Binomial Distributions: If the particle and anti-particle are independently distributed as Binomial or Negative Binomial distributions (BD/NBD). Then various order cumulants of the net-particle distributions can be expressed in term of cumulants of the particle and anti-particle distributions: $C_{n}^{net-p}=C_{n}^{p}+{(-1)^{n}}C_{n}^{\bar{p}}$. The first four order cumulants can be written as: $C_{2}^{x}=\sigma_{x}^{2}={\varepsilon_{x}}{\mu_{x}},C_{3}^{x}={S_{x}}\sigma_{x}^{3}={\varepsilon_{x}}{\mu_{x}}(2{\varepsilon_{x}}-1),C_{4}^{x}={\kappa_{x}}\sigma_{x}^{4}={\varepsilon_{x}}{\mu_{x}}(6\varepsilon_{x}^{2}-6{\varepsilon_{x}}+1)$ , where $\varepsilon_{x}=\sigma_{x}^{2}/\mu_{x}$, $\mu_{x}=M_{x}$, $M_{x}$ is the mean values of particles or anti-particles distributions, $x$=particle or anti-particle. $\varepsilon_{x}<1$ means the underlying distributions of particles or anti-particles are Binomial distributions, while $\varepsilon_{x}>1$ gives Negative Binomial distributions westfall . The input parameters for BD/NBD expectations are the measured mean and variance of the particle and anti-particle distributions.

IV.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 ShuHe_PLB . 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 ($\text{d}N/\text{d}y$) 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 BNS_2013_Koch ; HRG_FJH2 . 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 ($\Delta y$), the net-proton (baryon) cumulant ratios will decrease, reach a minimum and then increase, which is the typical effects of baryon number conservation BNS_2013_Koch ; HRG_FJH2 . 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 Kenji_MF .

IV.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 ampt and UrQMD urqmd 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 $\Delta$ resonance formation via $p\pi$ and $n\pi$ 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 Asakawa_formula . The converting formulas for various order net-baryon cumulants can be written as:

[TABLE]

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

IV.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 2015_Acceptance_Misha ; Corr_func_Koch_2016 , which are directly related to the correlation length ($\xi$) of system. The multi-particle density are related to factorial moments as:

[TABLE]

where $F_{n}$ is the $n^{th}$ order factorial moment and $\rho({p_{1}},...,{p_{n}})$ is the $n$ particle density distributions. The integral sums over the interested phase space. The generation function of the factorial cumulants is :

[TABLE]

[TABLE]

where $c_{k}$ is the $k^{th}$ order factorial cumulant, $N$ is the random variable. We have the relation between factorial moments and correlation function as:

[TABLE]

On the other hand, the relation between factorial cumulant ($c_{k}$) and cumulants can be expressed as:

[TABLE]

where the $s_{1}$ is the sterling number of the first kind, $C_{i}$ is $i^{th}$ order cumulant. Then, we have the following equations:

[TABLE]

[TABLE]

It is well known that high order cumulants ($C_{n},n>2$) are zero for the gaussian distribution and thus these are ideal probe of the Non-Gaussianity. For correlation function $c_{n}$ ($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 $\hat{c}_{k}$ as:

[TABLE]

where $k=2,3,4...,n$. The different order correlation strength parameter $\hat{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 :

[TABLE]

For eg., 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:

[TABLE]

The long range correlation become dominated near the critical point and the cumulant are dominated by the highest order correlation function as: $C_{k}\approx c_{k}\propto\langle N\rangle^{k}$. If the thermal statistical fluctuations dominated in the system, we have $C_{k}\approx c_{k}\propto\langle N\rangle$. 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, hadronic scattering.

Figure 15 shows the rapidity dependence of the proton cumulants and correlation functions in Au+Au collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}$=5 GeV from JAM transport model calculation JAM_model . 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.

V Data Analysis Methods

In the data analysis, we applied a series of analysis techniques to suppress backgrounds and make precise measurements of the event-by-event fluctuation analysis in heavy-ion collisions. Those include : (1) Centrality bin width correction WWND2011 ; technique . This is to remove centrality bin width effect, which is caused by volume variation within a finite centrality bin size. (2) Carefully define the collision centrality to suppress volume fluctuations and auto-correlations technique . (3) Efficiency correction for the cumulants. (4) Estimate the statistical error with Delta theorem and/or Bootstrap methods Delta_theory ; voker_eff1 ; voker_eff2 ; Unified_Errors . Those techniques are very crucial to precisely measure the dynamical fluctuation signals from heavy-ion collisions. Let’s discuss those techniques one by one.

V.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. Fig. 16 shows a Glauber Monte Carlo event of Au+Au collision at $\sqrt{s_{{}_{\mathrm{NN}}}}$ = 200 GeV with impact parameter $b=6$ fm miller2007glauber . The blue and red solid circles represent the participant nucleons from the two colliding gold nuclei. Fig. 17 shows the average number of participant nucleons ($\langle N_{part}\rangle$ ) and average number of binary collisions ($\langle N_{coll}\rangle$) 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  miller2007glauber . 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. Fig. 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.

V.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 WWND2011 ; technique . 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 ($C_{n}$) 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:

[TABLE]

where the $n_{r}$ is the number of events for multiplicity bin $r$ and the corresponding weight for the multiplicity $r$, ${\omega_{r}}={n_{r}}/\sum\limits_{r={N_{1}}}^{{N_{2}}}{{n_{r}}}$. $N_{1}$ and $N_{2}$ are the lowest and highest multiplicity values for one centrality bin. Once having the centrality bin width corrected cumulants via Eq.(86), we can calculate the various order cumulant ratios, for e.g. \kappa\sigma^{2}$$=C_{4}/C_{2} and S\sigma$$=C_{3}/C_{2}, where the $\kappa$ 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\sigma,\kappa\sigma^{2})$ of net-proton multiplicity distributions for Au+Au collisions at $\sqrt{s_{{}_{\mathrm{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 $n_{r}$ in Eq. (86) with $1/error^{2}$ 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 $\kappa\sigma^{2}$. 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.

V.3 Volume Fluctuations Effects

Volume fluctuations are long standing notorious background for the event-by-event fluctuation analysis in heavy-ion collisions technique ; Haojie_vol ; PBM_vol ; Skokov_vol ; haojie_plb ; bzdak2016correlated . 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 bzdak2016correlated 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_{{}_{\mathrm{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 $\sigma^{2}/M$ of $N_{part}$ distributions can be regarded as the centrality resolution for a certain centrality definition. Fig. 20 shows the centrality dependence of $\sigma^{2}/M$ of number of participant nucleons ($N_{part}$) distributions for Au+Au collisions at $\sqrt{s_{{}_{\mathrm{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 $\eta$ coverage ($|\eta|<0.5,1.0,1.5,2.0$). It shows that when we define the centrality with $|\eta|<2$, the $\sigma^{2}/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. Fig. 21 shows the energy dependence of moment product ($S\sigma,\kappa\sigma^{2}$) 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 $\eta$ 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 $\eta$ coverage of the charged particles used in centrality definition, we find the volume fluctuations are strongly suppressed. The $\kappa\sigma^{2}$ (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.

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. Fig.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 $\eta$ coverage ($|\eta|<$ 0.5 and 2). The larger $\eta$ coverage means more particles are included in the centrality definition and with better centrality resolution. Indeed, it shows that the results from wider $\eta$ coverage centrality definition are get suppressed comparing with the results from narrower $\eta$ coverage centrality definition. At fixed centrality definition, the results with CBWC are always smaller than the results without CBWC.

V.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 eg., 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. Fig. 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.

V.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 voker_eff1 ; voker_eff2 . Here, we provide a unified description of efficiency correction and error estimation for cumulants of multiplicity distributions Unified_Errors . 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 spectra_STAR . 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_{\bar{p}}$, which is proton number minus anti-proton number. The average value over the whole event ensemble is denoted by $\langle n\rangle$, 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 $n_{p}$ and anti-proton number $n_{\bar{p}}$ can be expressed as voker_eff1 :

[TABLE]

where the $P({N_{p}},{N_{\bar{p}}})$ is the original joint probability distribution of number of proton ($N_{p}$) and anti-proton ($N_{\bar{p}}$), $\varepsilon_{p}$ and $\varepsilon_{\bar{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:

[TABLE]

With the Eq. (87), (88) and (89), one can obtain a useful relation between the efficiency corrected and uncorrected factorial moments as:

[TABLE]

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 $\boldsymbol{X}=(X_{1},X_{2},...,X_{k})^{{}^{\prime}}$ and a set of number $\boldsymbol{r}=(r_{1},r_{2},...,r_{k})^{{}^{\prime}}$. The multivariate moments, central moments and factorial moments can be written as:

[TABLE]

where $E$ denotes the expectation value operator, and the ${m_{\boldsymbol{r}}}(\boldsymbol{X})$, ${\mu_{\boldsymbol{r}}}(\boldsymbol{X})$ and ${F_{\boldsymbol{r}}}(\boldsymbol{X})$ are multivariate moments, central moments and factorial moments, respectively. Then, we have the relation between the moments and central moments by using binomial expansions:

[TABLE]

where $\boldsymbol{i}=(i_{1},i_{2},...,i_{k})^{{}^{\prime}}$. To get the relation between moments and factorial moments, one needs the Stirling numbers of the first ($s_{1}(n,i)$) and second kind ($s_{2}(n,i)$), which are defined as:

[TABLE]

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

[TABLE]

and

[TABLE]

The Stirling number of the first kind may have the negative value while the value of the second kind is always non-negative. With the two kinds of Stirling numbers, one can write down the relations between moments and factorial moments as:

[TABLE]

[TABLE]

With Eq. (91) to (100), one can express the moments of net-proton distributions in terms of the factorial moments. There are two variables in net-proton number calculation, the number of protons ($N_{p}$) and anti-protons ($N_{\bar{p}}$). The $n^{th}$ order moments of net-proton distributions can be expressed in term of factorial moments:

[TABLE]

Actually, two steps are needed to obtain this equation, the first step is to expand the moments of net-proton to the bivariate moments by using binomial expansion, and the other one is to express the bivariate moments in term of the factorial moments using the Eq. (99). Now, one can easily calculate the efficiency corrected moments of net-proton distributions in heavy-ion collisions by using the Eq. (90) and (101). Finally, we can express the efficiency corrected cumulants of net-proton distribution with the efficiency corrected moments by using the recursion relation:

[TABLE]

where the $C_{r}$ denotes the $r^{th}$ order cumulants of net-proton distributions. In principle, one can also express the factorial moments in Eq. (101) in terms of the cumulants and the various order efficiency corrected cumulants can be expressed by the measured cumulants and efficiency as :

[TABLE]

where the $(X,Y)$ and $(x,y)$ are the numbers of $(p,\bar{p})$ produced and measured, respectively. $\varepsilon=\varepsilon_{p}=\varepsilon_{\bar{p}}$ is the $p(\bar{p})$ efficiency. Obviously, the efficiency corrected cumulants are sensitive to the efficiency and depend on the lower order measured cumulants. For more detail discussion of this method, one can also refer to Kitazawa_Eff .

In the previous discussion, the detection efficiency of proton and anti-proton are considered to be constant within the entire phase space. In many cases, the efficiency of proton and anti-proton will depend on the phase space (transverse momentum ($p_{T}$), rapidity (y), azimuthal angle ($\phi$)). In this sense, one has to re-consider the efficiency correction method. In the paper voker_eff2 , a new method for dealing with this case has been discussed, but the formulae for efficiency correction are rather involved and difficult to understand. In the following, we will provide an alternative efficiency correction method for the phase space dependent efficiency, which is straightforward and easier to understand. For simplify, we only consider the phase space of the proton and anti-proton are decomposed into two sub-phase spaces (1 and 2), within which the efficiency of proton and anti-proton are constant. We use the symbol $\varepsilon_{{p_{1}}},\varepsilon_{{p_{2}}}$ and $\varepsilon_{{{\bar{p}}_{1}}},\varepsilon_{{{\bar{p}}_{2}}}$ to denote the efficiency of proton and anti-proton in the two sub-phase spaces, and the corresponding number of proton and anti-proton in the two sub-phase spaces are $N_{p_{1}}$, $N_{p_{2}}$ and $N_{\bar{p}_{1}}$, $N_{\bar{p}_{2}}$, respectively. Using the relations in Eq. (99) and (100), one has:

[TABLE]

Based on the Eq. (104), we build up a relation between the bivariate factorial moments of proton and anti-proton distributions in the entire phase space and the multivariate factorial moments of proton and anti-proton distributions in the two sub-phase spaces. As a direct extension of Eq. (90) for multivariate case, the efficiency corrected multivariate factorial moments of proton and anti-proton distributions in the sub-phase spaces can be obtained as:

[TABLE]

where ${{f_{u,v,j,k}}(N_{{p_{1}}},N_{{p_{2}}},N_{{{\bar{p}}_{1}}},N_{{{\bar{p}}_{2}}})}$ is the measured multivariate factorial moments of proton and anti-proton distributions. By using Eq. (101), (102), (104) and (105), one can obtain the efficiency corrected moments of net-proton distributions for the case, where the proton (anti-proton) are with different efficiency in two sub-phase spaces. If the efficiency of the proton (anti-proton) have large variations within the phase space, one needs to further divide the phase space into small ones. It is easy and straightforward to do this, but it is time consuming and requires more computing resources.

To verify the phase space dependent efficiency correction formulas, we perform a calculation of the net-proton fluctuations with AMPT string melting model. The invariant $p_{T}$ spectra of proton and anti-proton from AMPT can be found in the Fig.24 left. In Fig.24 right, we set by hand the $p_{T}$ dependent efficiency for (anti-)protons with the efficiency at low $p_{T}$ ($0.4<p_{T}<0.8$ GeV/c) : 80% and high $p_{T}$ ($0.8<p_{T}<2$ GeV/c): 50%, respectively. The efficiency response function is set to be binomial distribution. Then, the measured net-proton distributions are the convolution between original model inputs and the binomial distributions. By doing this, we can calculate the measured cumulants of net-proton of Au+Au collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}$ = 39 GeV from AMPT string melting model with $p_{T}$ dependent efficiency. With the same procedures as we did in the real data analysis, we apply the phase space dependent efficiency formulas to do the efficiency correction for the measured cumulants. Fig. 25 shows that the efficiency corrected cumulants are consistent with the results from original model input within uncertainties. The statistical errors for the efficiency corrected cumulants are calculated from Delta theorem, which will be discussed later. Finally, this test confirms that the phase space dependent efficiency correction formulas we obtained are reliable and work well. On the other hand, if the efficiency response function is non-binomial type, instead of using analytical formulas, the unfolding method with real response matrix should be used Koch_eff_mul .

V.6 Error Estimation for the Efficiency Corrected Cumulants

Based on the Delta theorem in statistics, we obtained the error formulas for various order cumulants and cumulant ratios Delta_theory . However, those formulas can only be applied to the case, where the efficiency is unity ($\varepsilon=1$). It is not straightforward and easy to calculated the statistical errors for efficiency corrected cumulants with $\varepsilon\neq 1$ and one can not directly use the formulas obtained in the paper Delta_theory . In the following, we will derive general error formulas for estimating the statistical errors of efficiency corrected cumulants of conserved quantities in heavy-ion collisions based on the Delta theorem in statistics. With those analytical formulas, one can predict the expected errors with the number of events and efficiency numbers.

The Delta theorem in statistics is a fundamental theorem which is used to approximate the distribution of a transformation of a statistic in large samples if we can approximate the distribution of the statistic itself. Distributions of transformations of a statistic are of great importance in applications. We will give the theorem without proofs and one can see asytheory ; junshan .

Delta Theorem**:** Suppose that ${\bf{X}}=\{X_{1},X_{2},...,X_{k}\}$ is normally distributed as $N({\bf{\mu}},{\bf{\Sigma}}/n)$, with $\bf\Sigma$ a covariance matrix. Let ${\bf{g(x)}}=(g_{1}({\bf x}),...,g_{m}({\bf x}))$, ${\bf{x}}=(x_{1},...x_{k})$, be a vector-valued function for which each component function $g_{i}({\bf x})$ is real-valued and has a non-zero differential $g_{i}(\mu)$, at ${\bf{x}}={\bf{\mu}}$. Put

[TABLE]

Then

[TABLE]

where $n$ is the number of events.

Based on the Delta theorem, one can derive the general error formula for a statistic quantity. Suppose, statistic quantity $\phi$ is as a function of random variables ${\bf{X}}=\{X_{1},X_{2},...,X_{m}\}$, then the transformation functions ${\bf{g(X)}}=\phi({\bf X})$. The D matrix can be written as:

[TABLE]

and the covariance matrix $\Sigma$ is:

[TABLE]

Based on Eq. (107), the variance of the statistic $\phi$ can be calculated as:

[TABLE]

where $V(X_{i})$ is the variance of variable $X_{i}$ and $Cov(X_{i},X_{j})$ is the covariance between $X_{i}$ and $X_{j}$. To calculate the statistical errors, one needs to know the variance and covariance of the variable $X_{i}$ and $X_{j}$ in the Eq. (110). Since the efficiency corrected moments are expressed in terms of the factorial moments, the factorial moments are the random variable $X_{i}$ in Eq. (110). Then, we need to know the expression for variance and covariance of the factorial moments. It is known that the covariance of the multivariate moments advancetheory can be written as:

[TABLE]

where $n$ is the number of events, $m_{r,s}=<X_{1}^{r}X_{2}^{s}>$ and ${m_{u,v}}=<X_{1}^{u}X_{2}^{v}>$ are the multivariate moments, the $X_{1}$ and $X_{2}$ are random variables. Then, we can obtain the variance of the cumulants and cumulant ratios as:

[TABLE]

where $\mu_{r}=\langle(\delta N)^{r}\rangle$ is the $r^{th}$ order central moments, $m_{r}=\mu_{r}/\sigma^{r}$ and $n$ is the number of events. For normal distributions with width $\sigma$, the statistical error of the cumulants and cumulant ratios at different orders can be approximated as:

[TABLE]

Figure 26 shows the relative errors of cumulants and cumulant ratios of Skllellam distribution as a function of number of events $N$. It is found that the higher orders cumulants are with larger relative errors than the low orders at the same number of events $N$.

Based on Eq. (100) and (111), one can obtain the covariance for the multivariate factorial moments as:

[TABLE]

where the $f_{(r,u),(s,v)}$ is defined as:

[TABLE]

The definition of bivariate factorial moments $f_{r,s}$, $f_{u,v}$ and $f_{\alpha,\beta}$ are the same as Eq. (89). The Eq. (122) can be put into the standard error propagation formulae (110) to calculate the statistical errors of the efficiency corrected moments.

Besides the Delta theorem for estimating the statistical errors, another computer intensive one is the so called bootstrap, which is based on resampling method. with the bootstrap method, one needs to prepare $B$ new samples. Every new sample is sampling randomly with replacement from the original sample and are with the same number of events as the original one. The uncertainty on a statistic quantity is estimated by the root mean square of the $B$ values of the statistic quantity obtained from these samples. In the MC simulation, we set the number of new samples $B=200$. The variance of the statistic quantity $\Phi$ can be given by

[TABLE]

For comparison, we show the error estimation for the efficiency corrected $\kappa\sigma^{2}$ of Skellam distributions with Delta theorem and Bootstrap method in the Fig. 27 and Fig. 28, respectively. Both the Delta theorem and Bootstrap method can reasonably describe the statistical errors of the efficiency corrected $\kappa\sigma^{2}$ with various efficiency numbers ranging from 30% to 100%. The probability for the error bars of those data points touching the mean value are very close to the expected value 68%. Since we concentrate on the comparison of the magnitude of the statistical error calculated from the Delta theorem and Bootstrap methods, the data points are calculated from the same data sets and thus the $\kappa\sigma^{2}$ values are identical, while the statistical error bars of the data points in the two figures are not identical. This consistency verifies that the analytical error formulas derived from Delta theorem is correct. However, the calculation speed of Delta theorem method is much faster than that of Bootstrap method. On the other hand, since one cannot obtain events further into the tails than those in the original sample, the bootstrap method might run into difficulties if the quantity whose variance is being estimated depends heavily on the tails of distributions.

Figure 29 shows the statistical errors for the efficiency corrected $\kappa\sigma^{2}$, $S\sigma$ and $\sigma^{2}/M$ as a function of efficiency. In simulation, the efficiency effects are implemented for the original skellam distribution and the number of events is fixed to be one million for each data point. It can be found that the statistical errors are dramatically increase when decreasing the efficiency number, especially for higher order cumulant ratios. We also fit those data points with the functional form:

[TABLE]

where $n$ is the number of events which is fixed to be one million here, $a$ and $b$ are free parameters. The fitting results of $a$ and $b$ are 40.6 and 2.06 for $\kappa\sigma^{2}$, 6.02 and 1.65 for $S\sigma$, 4.96 and 0.89 for $\sigma^{2}/M$, respectively. The parameters $a$ and $b$ depend on the original distribution and the studied statistic quantity. We can understand the effects of the efficiency on the statistical errors in an intuitive way. The efficiency will cause the loss of information of the original distributions, especially at the tails. The smaller the efficiency is, larger uncertainties we will get for the efficiency corrected results and needs more events to recover the original information.

VI Experimental Results

One of the main goals of the beam energy scan program at RHIC is to explore the phase structure of the hot dense nuclear matter created in the relativistic heavy-ion collisions, especially searching for the QCD critical point and mapping out the first order phase boundary. From the year of 2010 to 2014, RHIC has finished the first phase of BES program, in which two gold nuclei collide at $\sqrt{s_{{}_{\mathrm{NN}}}}$ = 7.7, 11.5, 14.5 (taken at 2014), 19.6, 27, 39, 62.4 and 200 GeV. The STAR experiment has published the energy dependence of cumulants (up to fourth order) of net-proton 2010_NetP_PRL ; STAR_BES_PRL and net-charge netcharge_PRL multiplicity distributions in Au+Au collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}$= 7.7, 11.5, 19.6, 27, 39, 62.4 and 200 GeV. For net-proton analysis, the protons and anti-protons are identified with ionization energy loss in the Time Projection Chamber (TPC) of the STAR detector within the transverse momentum range $0.4<p_{T}<0.8$ GeV/c and at mid-rapidity $|y|<0.5$. For the net-charge, the charged particles are measured within transverse momentum range $0.2<p_{T}<2$ GeV/c and pseudo-rapidity range $|\eta|<0.5$.

Figure 30 shows the energy dependence of cumulant ratios of net-proton and net-charge distributions of Au+Au collisions for two centralities (0-5% and 70%-80%) at $\sqrt{s_{{}_{\mathrm{NN}}}}$= 7.7, 11.5, 19.6, 27, 39, 62.4 and 200 GeV. The Skellam (Poisson) expectations shown in the figure reflect a system of totally uncorrelated, statistically random particle production. It predicts the $\kappa\sigma^{2}$ and $S\sigma$/Skellam to be unity for Skellam expectations as well as in the hadron resonance gas model. For the net-proton results, the most significant deviation of $S\sigma$ and $\kappa\sigma^{2}$ from Skellam distribution is observed at 19.6 and 27 GeV for 0-5% Au+Au collisions. At energies above 39 GeV, the results are close to Skellam expectation. As the statistical errors are large at low energies (7.7 and 11.5 GeV), more statistics is necessary to quantitatively understand the energy dependence of $S\sigma$ and $\kappa\sigma^{2}$. To understand the effects of baryon number conservation etc., UrQMD model calculations (a transport model which does not include a CP) for 0-5% are presented and the results show a monotonic decrease with decreasing beam energy. For more details on baseline comparison, one can see QM2014_baseline . For the net-charge results, we did not observe non-monotonic behavior for $S\sigma$ and $\kappa\sigma^{2}$ within current statistics. The expectations from negative binomial distribution can better describe the net-charge data than the Poisson (Skellam) distribution. More statistics is needed for the measurements of net-charge moments.

In the CPOD2014 2014_Luo_CPOD and QM2015 conferences QM_XFLUO ; 2015_Jochen_QM , the STAR experiment reported the preliminary results of net-proton fluctuations with wider transverse momentum coverage ($0.4<p_{T}<2$ GeV/c). In the new results, the $p_{T}$ range of (anti-)protons are extended from $0.4<p_{T}<0.8$ to $0.4<p_{T}<2$ GeV/c. This is realized by using the Time of Flight (ToF) detector to identify the high $p_{T}$ ($0.8<p_{T}<2$ GeV/c) (anti-)protons. At low $p_{T}$ region ($0.4<p_{T}<0.8$ GeV/c), only Time Projection Chamber (TPC) is used to identify the (anti-)protons whereas the (anti-)protons at high $p_{T}$ ($0.8<p_{T}<2$ GeV/c) are jointly identified by TPC and ToF. Fig. 31 left show the particle identification (PID) plot for TPC and ToF detector. The white dashed boxes in the ToF PID plot denote the protons (upper) and kaons (lower) PID cuts region, respectively. Fig. 31 right show the proton phase space in Au+Au collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}$ = 14.5 GeV measured by the STAR experiment. The protons and anti-protons in the regions covered by the blue dashed boxes are used in the net-proton fluctuation analysis. Figure 32 shows the uncorrected event-by-event net-charge, net-kaon and net-proton multiplicity distributions in three centralities (0-5%, 30-40% and 70-80%) for Au+Au collisions at $\sqrt{s_{NN}}$ = 14.5 GeV. Those raw distributions from a wide centrality bin can not be used to calculate the various order cumulant directly due to the effects of finite efficiency and volume variation. However, there are some theoretically works about using those distributions to extract criticality kenji ; kenji_netp ; PBM_netq ; PBM_netpdis . The shape of net-particle multiplicity distributions for different centralities are different. The standard deviation $\sigma$ of the net-particle distributions get bigger for central collisions than peripheral and mid-central. We also observed that the net-charge multiplicity distributions have the largest standard deviation, $\sigma$, comparing with the net-proton and net-kaon distributions at fixed centrality. As shown in Eq. (120), the statistical errors of the $r^{th}$ order cumulants are proportional to the $r^{th}$ power of the standard deviation ($\sigma^{r}$). This indicates that with the same number of events, the net-charge fluctuations measurements will have much lager statistical errors than the results of net-proton and net-kaon fluctuations. Detailed discussions about the efficiency correction and error estimation can be found in technique ; Unified_Errors .

Figure 33 shows the centrality dependence of detection efficiency for (anti-)protons in two $p_{T}$ ranges ( $0.4<p_{T}<0.8$ and $0.8<p_{T}<2$ GeV/c) in Au+Au collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}$=7.7 , 11.5, 19.6, 27, 39, 62.4 and 200 GeV. The efficiency of protons and anti-protons at high $p_{T}$, $0.8<p_{T}<2$ GeV/c is smaller than that of low $p_{T}$, $0.4<p_{T}<0.8$ GeV/c. This is because, besides the time projection chamber (TPC), the time of flight (ToF) detector is used to identify the high $p_{T}$ (anti-)protons and the ToF matching efficiency is introduced in addition to the TPC tracking/acceptance efficiencies. While at low $p_{T}$, only TPC is used to identify protons and anti-protons. Thus, the average efficiency for protons or anti-protons at low $p_{T}$ and high $p_{T}$ can be calculated as:

[TABLE]

where the $\varepsilon({p_{T}})={\varepsilon_{tpc}}({p_{T}})$ for $0.4<p_{T}<0.8$ GeV/c and $\varepsilon({p_{T}})={\varepsilon_{tpc}}({p_{T}}){\varepsilon_{tof}}({p_{T}})$ for $0.8<p_{T}<2$ GeV/c. The efficiency corrected $p_{T}$ distribution function $f({p_{T}})$ is defined as $f({p_{T}})=dN/d{p_{T}}$. The TPC efficiency ($\varepsilon_{tpc}({p_{T}})$) of protons or anti-protons are obtained from the so-called embedding simulation techniques and the ToF matching efficiency ($\varepsilon_{tof}({p_{T}})$) can be calculated from the real data. The average efficiencies of protons and anti-protons have centrality (multiplicity) dependence and increase from central to peripheral collisions for all energies. Due to material absorption of anti-protons in the detector, the efficiencies of anti-protons are always slightly lower than protons.

Figure 34 shows the centrality dependence of efficiency corrected cumulants ($C_{1}\sim C_{4}$) of net-proton, proton and anti-proton distributions in Au+Au collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}$=7.7, 11.5, 14.5, 19.6, 27, 39, 62.4 and 200 GeV. The protons and anti-protons are measured within transverse momentum $0.4<p_{T}<2$ GeV/c and at mid-rapidity ($|y|<0.5$). At high energies, the cumulants (up to fourth order) of net-proton, proton and anti-proton distributions show a linear dependence on the average number of participant nucleons ($\langle N_{part}\rangle$). This is consistent with the additive properties of the cumulants that the system consists of many multi-independent emission sources of protons and anti-protons and those emission sources are linear dependent on the system volume (centralities). The proton cumulants are always larger than the anti-proton cumulants and the difference between proton and anti-proton cumulants are larger in low energies than high energies. The cumulants of net-proton distributions closely follow the proton cumulants when the colliding energy decreases. These observations can be explained as the interplay between the baryon stopping and pair production of protons and anti-protons. At high energies, protons and anti-protons mainly come from the pair production and the number of protons and anti-protons are very similar. At low energies, the production of protons is dominated by baryon stopping and the number of protons is much larger than the number of anti-protons. The efficiency corrected fourth order net-proton and proton cumulants ($C_{4}$) of 7.7 and 11.5 GeV significantly increase in the $0\sim 5\%$ and $5\sim 10\%$ centrality bins with respect to the efficiency uncorrected results. It means the efficiency corrections are big effects, especially for the high order cumulants. Furthermore, the efficiency correction not only affects the values but also lead to increasing of the statistical errors for the various order cumulants, as $error(C_{n})\sim\sigma^{n}/\varepsilon^{\alpha}$, where the $\sigma$ in numerator is the standard deviation of the particle distributions and the denominator $\varepsilon$ is the efficiency number, $\alpha$ is a positive real number Unified_Errors .

The STAR Collaboration reported preliminary results of cumulants of net-kaon distributions in Au+Au collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}$= 7.7, 11.5, 14.5, 19.6, 27, 39, 62.4 and 200 GeV in the QM2015 conference 2015_Jochen_QM . The net-kaon fluctuations is used to approximate the fluctuations of net-strangeness, a conserved charge in strong interaction. The susceptibilities of net-strangeness can be computed in Lattice QCD. The $K^{+}$ and $K^{-}$ are measured with transverse momentum $0.2<p_{T}<1.6$ GeV/c and at mid-rapidity $|y|<0.5$. At low $p_{T}$ region ($0.2<p_{T}<0.4$ GeV/c), the charged kaons are identified by TPC only whereas at high $p_{T}$ ($0.4<p_{T}<2$ GeV/c), ToF is also used in addition with TPC. To avoid auto-correlation, the collision centrality is determined by measured charged particles within $|\eta|<1$ excluding charged kaons. Fig. 35 shows the efficiency corrected centrality dependence of cumulants ($C_{1}\sim C_{4}$) of net-kaon multiplicity distributions in Au+Au collisions at $\sqrt{s_{\rm NN}}$ =7.7$\sim$200 GeV. The red dashed lines represent the Poisson expectations, where the probability distributions of the $K^{+}$ and $K^{-}$ are assumed to be the independent Poisson distributions. In general, various order cumulants show a linear variation with the averaged number of participant nucleons ($\langle N_{part}\rangle$). The variance are systematically below the Poisson expectations, especially at high energies. It means that the $K^{+}$ and $K^{-}$ are correlated with each other due to the pair productions. However, the $C_{3}$ and $C_{4}$ are consistent with Poisson expectation within uncertainties. The large uncertainties observed in the $C_{3}$ and $C_{4}$ are due to the low detection efficiency of kaons ($\sim 40\%$).

Figure 36 left shows the energy dependence of cumulants ($C_{1}\sim C_{4}$) for net-kaon, $K^{+}$, and $K^{-}$ multiplicity distributions in Au+Au collisions measured by the STAR experiment. The mean values of the $K^{+}$ and $K^{-}$ show monotonic decreasing trends when the energy decrease. Furthermore, the mean values of $K^{+}$ is always above $K^{-}$, and the difference between these two values are bigger at lower energies. These two observations are due to interplay of the pair and associate production for $K^{+}$, and $K^{-}$ as a function of collisions energies. In addition to the pair production of $K^{+}$ and $K^{-}$, the $K^{+}$ is also produced by the associate production with $\Lambda$ hyperon and the fraction of $K^{+}$ from associate production is lager at low energies than at high energies. It also leads to the increasing of the net-kaon mean values when decreasing the energies. The corresponding Poisson expectations are also plotted as different lines for comparison. In general, the cumulants of $K^{+}$, and $K^{-}$ distributions are consistent with the Poisson baseline within uncertainties. Due to the correlation between $K^{+}$ and $K^{-}$ , the variance of the net-kaon distributions are smaller than its Poisson expectations, in which one assumes the independent of the $K^{+}$ and $K^{-}$. The higher order net-kaon cumulants are consistent with Poisson expectations within uncertainties. Fig. 36 right shows the energy dependence of cumulants ($C_{1}\sim C_{4}$) of net-proton, proton and anti-proton multiplicity distributions in Au+Au collisions measured by the STAR experiment. The mean values of protons and net-protons show monotonic increasing trends when decreasing the colliding energy whereas the mean values of anti-protons show opposite trend. Those can be understood in terms of the interplay between the baryon stopping and pair production for proton and anti-proton as a function of collision energy. At low energies, the baryon stopping becomes more dominate while at high energies, the pair production is the main production mechanism of the proton and anti-protons. In the figure, it also shows the comparison between the cumulants of net-proton, proton and anti-proton distributions and the corresponding Poisson expectations. We found that the higher the order of the cumulant, the larger the deviations from the Poisson expectation for the net-proton and proton. Largest deviations are found for $C_{4}$ at 7.7 GeV. The cumulants of anti-proton distributions can be described by the Poisson expectations very well. More baselines discussions from Hadronic Resonance Gas model, transport model UrQMD, binomial and negative binomial have been also discussed.

Figure 37 shows the energy dependence of cumulant ratios ($\sigma^{2}/M$, $S\sigma/$Skellam, $\kappa\sigma^{2}$) of net-charge WWND2016_Xuj ; 2015_Jochen_QM , net-kaon WWND2016_Xuj ; SQM2016_Xuj , and net-proton 2014_Luo_CPOD multiplicity distributions in Au+Au collisions measured by the STAR experiment. The black solid circles on the left figure represent the results from Au+Au collisions at $\sqrt{s_{NN}}$ = 14.5 GeV, which is taken in the year 2014 and added into the trend of the published net-charge results netcharge_PRL (open stars). The bands are the results from UrQMD calculations without including the critical physics. The Poisson expectations are displayed as dashed lines. The $S\sigma$ values have been normalized by the Poisson expectations, the Skellam distributions. Thus, the Poisson expectations for both the $S\sigma/$Skellam and $\kappa\sigma^{2}$ are unity. It can be found that the $\sigma^{2}/M$ of net-charge, net-kaon and net-proton monotonically increase when increasing the collision energy. On the other hand, both the $S\sigma/$Skellam and $\kappa\sigma^{2}$ show weak energy dependence for net-charge and net-kaon measurements. No significant deviations from the Poisson expectations and UrQMD calculation are observed for net-charge and net-kaon cumulant ratios $S\sigma/$Skellam and $\kappa\sigma^{2}$ within uncertainties. We observe a clear non-monotonic energy dependence of net-proton $\kappa\sigma^{2}$ in top 0-5% central Au+Au collisions. The 0-5% net-proton $\kappa\sigma^{2}$ values are close to unity for energies above 39 GeV and show large deviations below unity around 19.6 and 27 GeV, and then increasing above unity below 19.6 GeV. The UrQMD calculations of net-proton $\kappa\sigma^{2}$ displaying a strong suppression below unity at lower energies is due to the effects of baryon number conservation. However, this suppression are not observed at low energies in the STAR data.

Figure 38 panels (a), (c), (d) show the energy dependence of $\kappa\sigma^{2}$ of net-charge, net-kaon and net-proton multiplicity distributions in Au+Au collisions measured by the STAR experiment for two centralities (0-5% and 70%-80%) at $\sqrt{s_{{}_{\mathrm{NN}}}}$= 7.7, 11.5, 14.5, 19.6, 27, 39, 62.4 and 200 GeV. The $\kappa\sigma^{2}$ of net-charge distributions in Au+Au collisions $\sqrt{s_{{}_{\mathrm{NN}}}}$= 7.7, 19.6, 27, 39, 62.4 and 200 GeV measured by the PHENIX experiment NetQ_PHENIX are shown in the panel (b) of Fig. 38. It is observed that the $\kappa\sigma^{2}$ of the net-charge and net-kaon distributions measured by the STAR experiment are with larger statistical errors than the errors of net-proton $\kappa\sigma^{2}$. This is because the statistical errors of $\kappa\sigma^{2}$ depend on the width ($\sigma$) of the multiplicity distributions (error({{\mbox{\kappa\sigma^{2}}}})\propto\sigma^{n-2}/(\sqrt{N}\epsilon^{n})) and the net-charge distributions are much wider than the net-proton and net-kaon. On the other hand, due to decay of kaons, the efficiency of kaon ($\sim 40\%$) is much lower than proton ($\sim 80\%$), this also leads to larger statistical errors for net-kaon fluctuations. For the STAR net-charge and net-kaon results, we did not observe non-monotonic behavior for $\kappa\sigma^{2}$ within current statistics. The Poisson expectations shown as dashed lines in the figure with unity value reflect a system of totally uncorrelated, statistically random particle production. It predicts the $\kappa\sigma^{2}$ to be unity for Poisson expectations as well as in the hadron resonance gas model. However, the expectations from negative binomial distribution can better describe the net-charge and net-kaon data than the Poisson expectations. The PHENIX net-charge $\kappa\sigma^{2}$ are with smaller errors than the results measured by the STAR experiment. This is because the PHENIX detector has much smaller acceptance than the STAR detector and thus the width of the net-charge distributions measured by the PHENIX experiment is much narrower. We observe a clear non-monotonic energy dependence for net-proton $\kappa\sigma^{2}$ in the most central ( 0-5%) Au+Au collisions with a minimum around 19.6 and 27 GeV. This non-monotonic behavior cannot be described by various model calculations without including CP physics Xuji_PRC ; ShuHe_PLB . Another model calculation with volume fluctuations also failed to describe this increasing at low energies bzdak2016correlated . At energies above 39 GeV, the 0-5% net-proton $\kappa\sigma^{2}$ are close to Poisson expectations while at energies below 19.6 GeV, it shows large increasing above unity. This large increase in 0-5% net-proton $\kappa\sigma^{2}$ at low energies

We want to make several remarks : (1) One needs to remember that the resonance decay effects are not excluded in the current experimental measurements of fluctuations of net-proton, net-kaon and net-charge. Based on the hadron resonance gas model calculation HRG_baseline , the decay effects for net-proton $\kappa\sigma^{2}$ is small and at 2% level. While for the net-charge, the decay effects are large. (2) The statistical error of cumulants ($\Delta(C_{n}$)) are related to the width of the distribution as $\Delta(C_{n})$ $\sim$ O($\sigma^{n}$) Delta_theory ; Unified_Errors . Thus, the wider is the distribution, the larger are statistical errors for the same number of events. (4) It is predicted by theoretical calculations that the net-baryon fluctuations are more sensitive to the criticality than the net-charge and net-strangeness Wenkai_NJL ; ratioCumulant . (5) The measurements of fluctuations of conserved quantities can be used to determine the freeze-out conditions in heavy-ion collisions by comparing with the Lattice calculations and/or HRG calculations bazavov2012freeze ; Borsanyi_freeze ; KaonFlu_Lattice ; Alba_freeze .

VII Beam Energy Scan Phase-II and STAR Detector upgrades

A second phase of the beam energy scan (BES-II) program at RHIC has been planned in the years 2019-2020 and focusing on energy rang $\sqrt{s_{{}_{\mathrm{NN}}}}$=$7.7\sim 20$ GeV BESII_WhitePaper . The long beam bunches and stochastic electron cooling technique will be used to accelerate gold beams, which will increase the luminosity about by a factor of 5-15 for corresponding collision energies compared to the BES-I. Since the luminosity will decrease as decreasing the colliding energy, the increasing of the luminosity is much more necessary and important at low energies, such as 7.7 GeV. This enable us to collect more number of events ($\sim 10-20$ times) to confirm the non-monotonic trends observed in the BES-I data. Furthermore, it will allow us to draw a solid conclusion and have more complete physical pictures from various experimental measurements. To study the QCD phase structure at high baryon density, operating the STAR detector at fixed-target mode has been also proposed. In the BES-II, fixed-target mode Au+Au collisions allow us to have energy coverage from $\sqrt{s_{{}_{\mathrm{NN}}}}$=3 GeV ($\mu_{B}$=720 MeV) up to 7.7 GeV. In Fig.39 left, the inner TPC (iTPC) of STAR is to be upgraded to improve the energy loss resolution and can extend the pseudo-rapidity coverage from $|\eta|<1$ to $|\eta|<1.5$ iTPC_proposal . It is also planed to install an end cap Time-of-Flight (eTOF) detector at the west end of the STAR TPC to extend the PID capability in the forward region eTOF_STAR . The iTPC upgrade is very important to search for the criticality and study the dynamical evolution of the fluctuations by looking at the rapidity acceptance dependence for the fluctuations of the conserved quantities 2015_Kitazawa_Rapidity ; Asakawa_review . In the forward and backward region of STAR detector, a new Event Plane Detector (EPD) will be also built and used to replace the Beam-Beam Counters (BBC) detector for centrality and event plane determination, which can be used to suppress the volume fluctuation and auto-correlation in the fluctuation analysis. In Fig.39 right, the blue band is the extrapolating from current measurements by assuming a power law behavior induced by critical fluctuations ({\mbox{\kappa\sigma^{2}}}{\propto}N^{3} 2015_Acceptance_Misha ). The width of the blue band is the estimated statistical errors with the BES-II statistics and iTPC upgrades.

Figure 40 shows the STAR preliminary results of energy dependence of the fourth-order fluctuations ($\kappa\sigma^{2}$) of net-proton, proton and anti-proton from the most top 5% central Au+Au collisions. Those data were taken from the first phase of the RHIC beam energy scan (BES-I) and from the kinematic region of mid-rapidity $|y|<0.5$ and transverse momentum $0.4<p_{T}<2$ GeV/c. Non-monotonic energy dependence is clearly shown in the $\kappa\sigma^{2}$ of net-proton and proton distributions. Although the statistical errors are large, the data shows a strong enhancement at the highest $\mu_{B}\sim 420$ MeV, corresponding to the Au+Au central collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}$ = 7.7 GeV. This indicates an attractive correlation in nature at the large baryon density region. On the other hand, the results from the transport model UrQMD (yellow-line) show a monotonic decrease from low to high baryon density region, reflecting the fact that the baryon number conservation in such high-energy nuclear collisions. All known model calculations have shown just that. It appears that the baryon number conservation is dominant in those model simulations. Note that in the Poisson limit, the absence of criticality or other dynamical correlations, the $\kappa\sigma^{2}$ is expected to be unity.

The green region in the figure is the projected error of the fourth order fluctuations $\kappa\sigma^{2}$ of net-protons in the second phase of the RHIC Beam Energy Scan (BES-II) program BESII_WhitePaper . The BES-II program, which is scheduled to take place during the years of 2019 and 2020 for the Au+Au collisions at 7.7-19.6 GeV, will take about 10 to 20 times (depending on energy) higher statistics data to confirm the non-monotonic behavior observed in the fourth order fluctuations ($\kappa\sigma^{2}$) of net-protons in Au+Au collisions in the BES-I measured by STAR. Since no one expects protons freeze-out at the critical point so experimentally one should search for the critical region instead of a point ratioCumulant ; Neg_Kurtosis . Assuming the data in the figure is related to the critical region, one must study the net-proton fluctuations at even higher baryon density region, i.e. at lower collision energies. At energy below 7.7 GeV, the collider mode experiments become inefficient so the fixed-target (FXT) mode is the way out.

VIII Summary

In this review, we summarized the fluctuations (up to fourth order) of net-proton, net-charge and net-kaon in Au+Au collisions at $\sqrt{s_{{}_{\mathrm{NN}}}}$= 7.7, 11.5, 14.5, 19.6, 27, 39, 62.4 and 200 GeV. Those data are taken in the year 2010 to 2014 and in the first phase beam energy scan program at RHIC. The corresponding baryon chemical potential ($\mu_{B}$) coverage is from about 23$\sim$420 MeV. To make precise measurements, a series of data analysis techniques have been built up to suppress the volume fluctuation and auto-correlation backgrounds. We also provide a unified description of the finite detection efficiency correction and error estimation for the various order cumulants of net-particle distributions. The statistical errors of the cumulants are related to the measured standard deviation of distributions ($\sigma$) and the particle detection efficiency as $error({C_{n}})\propto\sigma^{n}/(\sqrt{N}\epsilon^{n})$.

In summary, we have:

Experimental Observations:

1. (1)

Due to larger width of the net-charge distribution and lower efficiency of charged kaons, we have bigger statistical errors of cumulants of net-charge and net-kaon than the net-proton cumulants. Within current statistical uncertainties, the energy dependence of the net-charge and net-kaon $S\sigma$ and $\kappa\sigma^{2}$ are flat and consistent with Poisson expectations and UrQMD model calculations. 2. (2)

In general, various order cumulants show linear variation with the average number of participant nucleons ($\langle N_{part}\rangle$). The interplay of the production mechanisms for particle and anti-particle as a function of collision energy have significant impacts on the energy dependence of the cumulants. 3. (3)

We observed a clear non-monotonic energy dependence for the $\kappa\sigma^{2}$ of the net-proton, proton multiplicity distributions in 0-5% most central Au+Au collisions measured by the STAR experiment.

Theoretical and Model Calculations:

1. (1)

The non-monotonic behavior observed in the energy dependence of the 0-5% net-proton, proton $\kappa\sigma^{2}$ in Au+Au collisions are consistent with a presence of QCD critical point from model calculations, such as $\sigma$ field, NJL , PNJL and PQM models. Those model calculations suggested an non-monotonic oscillation pattern due to the sign change of the critical contributions in different QCD critical regions. However, it is still not conclusive yet and more works about the dynamical modelling of heavy-ion collisions are needed. 2. (2)

The large increasing in the net-proton and proton $\kappa\sigma^{2}$ at low energies cannot be reproduced by various transport model calculations. All known transport model calculations show a strong suppression with respect to unity at low energies, which is dominated by the effects of baryon number conservations.

Future directions:

1. (1)

Experimentally, in order to confirm the observed energy dependence structures in the high moments of net-protons in BES-I, the second phase of the beam energy scan (BES-II) at RHIC has been planned in 2019-2020 with increased luminosity BESII_WhitePaper . This allows us to have 10 to 20 times more statistics at energies $\sqrt{s_{{}_{\mathrm{NN}}}}$ = 7.7$\sim$19.6 GeV to explore the phase structure this low energy range with high precision. The upgrades of iTPC and EPD are ongoing in the STAR and will provide large rapidity coverage and forward centrality determination in the BES-II, respectively. The large rapidity coverage is very important for us to perform the rapidity dependence for the fluctuation analysis, which is crucial to test the long range correlation as well as power law behavior induced by QCD critical point. For energies below 7.7 GeV, the fixed-target mode becomes more efficient than collider mode. It also has been tested and proposed to operate the STAR detector under a fixed target mode in the BES-II. A fixed-target experiment called Compressed Baryonic Experiment (CBM) at FAIR CBM_physics will start in 2024 and, in its first phase SIS100, will cover the Au+Au collision energy range of $\sqrt{s_{{}_{\mathrm{NN}}}}$ = $2.5\sim 4.7$ GeV. This will be an ideal experiment to search for the QCD critical point at the high-baryon density region with high precision. 2. (2)

We have mentioned that the first-order phase boundary, the critical point and the smooth crossover are closely related thermodynamically. At high net-baryon density region, we are searching for the signatures of the QCD critical point and/or the first-order phase boundary. However, in the near future, at the high-energy frontier, one should also search for the experimental evidence of the smooth crossover. This can be done with higher order fluctuations of conserved quantities. At the vanishing baryon chemical potential, $\mu_{B}\sim 0$, although the transition is a smooth crossover, there should have the remnant criticality of the chiral transition. Higher order fluctuations, cumulants $C_{6}$ (sixth order) or $C_{8}$ (eighth order) could show strong oscillation and should be able to pick up the possible signal in heavy-ion collisions at both RHIC and LHC. These results will not only confirm experimentally the smooth crossover nature of the transition, may also provide the information on the width of the crossover, which is one of the key information of the QCD phase diagram at small net-baryon density. On the other hand, the measurements of the various order correlation functions as a function of centrality, rapidity and energy are also very useful to further understand the critical and non-critical physics contributions. 3. (3)

We also want to point out that due to density fluctuations near the QCD critical point, light nuclei production and/or nucleon-clusters, such as deuteron, 3He and 4He as well as the energy dependence of the low mass di-lepton yield CBM_physics ; rapp2013dilepton ; Yifei_dilepton ; rapp2002chiral could also be used to aid and complement to the critical point searches at the high baryon density region. Of course these different observables are with different systematics. Details analysis are needed in order to understand these systematic effects. 4. (4)

Theoretically, careful modellings for the critical fluctuations and dynamical evolution of the thermodynamic medium created in the heavy-ion collision at different energies are needed to understand the phase structure of QCD, in particular the de-confinement transition and possible critical point. Many attempts and progress have been made by physicist worldwide 2015_Swagato_evolution ; Swagato_2016_PRL ; Huichao_flu_QM2015 ; 2016_fluctuations_Huichao . Those theoretical inputs are particularly important to establish definitive connections between experimental observables and phase structures in the QCD phase diagram.

Acknowledgement

The work was supported in part by the MoST of China 973-Project No.2015CB856901, NSFC under grant No. 11575069.