{'fraction_non_alphanumeric': 0.033790745384089355, 'fraction_numerical': 0.013995896968315477, 'mean_word_length': 4.719440967876512, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Search for the QCD Critical Point in High Energy Nuclear Collisions

QCD critical point is a landmark region in the QCD phase diagram outlined by temperature as a function of baryon chemical potential. To the right of this second-order phase transition point, one expects first order quark-hadron phase transition boundary, towards the left a crossover region, top of it lies the quark gluon plasma phase and below it the hadronic phase. Hence locating the QCD critical point through relativistic heavy-ion collision experiments is an active area of research. Cumulants of conserved quantities in strong interaction, such as net-baryon, net-charge, and net-strangeness, are suggested to be sensitive to the physics of QCD critical point and are therefore useful observables in the study of the phase transition between quark-gluon plasma and hadronic matter. We review the experimental status of the search for the QCD critical point via the measurements of cumulants of net-particle distributions in heavy ion collisions. We discuss various experimental challenges and associated corrections in such fluctuation measurements. We also comment on the physics implications of the measurements by comparing them with theoretical calculations. This is followed by a discussion on future experiments and measurements related to high baryonic density QCD matter.


Introduction: The phase diagram of the QCD and heavyion collisions
Quantum Chromodynamics (QCD) is the underlying theory of strong interaction where the fundamental constituents are quarks and gluons. Under high temperature and pressure, the nuclear matter undergoes a transition from a confined state of quarks and gluons called the hadronic phase [1,2,3], to a state where quarks and gluons are deconfined, known as the Quark-Gluon Plasma (QGP) [4]. Experimental study of heavy-ion collisions has revealed the strongly coupled nature of this state of matter [5,6,7]. The investigation of the phase diagram of QCD is one of the primary goals of the heavy-ion collision experiments and has been an active subject of research for many years. Unlike the phase diagram of water which has gathered a competent understanding having electromagnetism as the nature of underlying interactions, the QCD phase diagram remains largely conjectured. QCD phase diagram can be represented by a graph showing the variation of temperature (T ) vs. chemical potential (µ) associated with conserved charges like baryon number (B), electric charge (Q), and strangeness number (S). In heavy-ion collisions, the value of µ Q and µ S are fixed by the ratio of electric charge to baryon number and net-strangeness content of the colliding nuclei. Thus, the experimentally accessible QCD phase diagram effectively reduces to a two-dimensional graph of T vs. µ B . A conjectured QCD phase diagram is shown in Fig. 1.
QCD calculations on lattice at vanishing baryonic chemical potential and large temperature have established that the nature of quark-hadron phase transition is a smooth crossover [8]. In lattice-QCD calculations, the chiral susceptibility is the order parameter distinguishing the phases. It is defined as, χ(N s , N τ ) = (∂ 2 /∂m 2 ud )(T /V ) log Z where Z is the partition function, m ud is the mass of the light u,d quarks, N s is the spatial extension, N τ is euclidean time extension, and V the system volume in lattice-QCD. The smooth crossover was proved by showing that the temperature dependence of the peak and the width of the chiral susceptibility are independent of the system volume. For a typical first-order phase transition, the height of the susceptibility peak should have been proportional to the volume, and the width would vary inversely with volume. For a second-order transition, a singular behavior should have been observed with the volume of the system (V α , α is a critical exponent). The point of sharpest change in the temperature dependence of the chiral susceptibility, the strange quark number susceptibility, and the renormalized Polyakov-loop (an approximate order parameter for quark deconfinement in a hot gluonic medium [9] defined as L(T ) ∼ lim r→∞ exp (−V (r)/T ), where V (r) is the potential between a static quark-antiquark pair separated by a distance r) are used to estimate the QCD transition temperature in lattice calculations. The crossover transition temperature at µ B = 0 is reported to be 156.5 ± 1.5 MeV [10]. The chiral phase transition temperature at µ B = 0, which is expected to put an upper bound on the temperature for the possible critical point at nonzero baryon chemical potential, has been recently reported to be 136 +3 −6 MeV [11] for vanishing two light quark masses and with the strange quark mass fixed to its physical value. The temperature reduces to 98 +3
Various QCD-based model calculations at large µ B predict the phase transition to be of first order [13,14]. The point on the T vs. µ B graph, where the first-order phase transition changes to crossover, is called the QCD critical point (CP) [15,16,17]. The existence of the QCD critical point guarantees the occurrence of the first-order phase transition, and the vice-versa also holds true. There are several theoretical challenges in determining the location of the QCD critical point. The first principle lattice calculations at finite µ B faces the issue of sign-problem [18]. Even though the location of QCD critical point remains highly uncertain, recent estimates from the lattice-QCD calculations suggest, critical point, if it exists, should be in the region µ B ≥ 300 MeV [19].  : Chemical freeze-out parameters T and µ B obtained by thermal model fits to experimental data on particle yields as a function of collision energies [21]. Red line represents parametrisations of T and µ B as a function of √ s N N taken from Ref. [20] Experimentally, changing the energy and species of colliding nuclei can probe various regions of the QCD phase diagram. The measurements on particle yields and their ratios in the heavy-ion collision experiments can be compared to statistical thermal models to extract T and µ B at the chemical freezeout [20,21,22]. Figure 2 shows the extracted freeze-out parameters as a function of collision energy. The measurements suggest that µ B decreases with increasing center-of-mass energy ( √ s N N ) while T exhibits the opposite dependence for the lower collision energies ( √ s N N 20 GeV) then saturates towards higher √ s N N . Thus by varying the collision energy of heavy-ions the T and µ B axes of the phase diagram can be scanned. Varying the impact parameter of the collisions, the species of colliding ions, and the rapidity acceptance could possibly further broaden the experimentally accessible region of the phase diagram.
In the scenario where freeze-out happens close to the phase transition, event-by-event fluctuations of conserved charge multiplicity distributions like the net-baryon, net-charge, and net-strangeness are suggested as sensitive observables to explore the QCD phase structure and search for the QCD critical point [16,23]. The cumulants of event-by-event conserved charge multiplicity distributions are related to the correlation length (ξ) of the hot and dense medium formed in heavy-ion collision experiments. Upon approaching the critical point, the correlation length should diverge for a static and infinite system. This is also true for any response functions like the susceptibilities, compressibilities or heat capacities. However, the matter formed in heavy-ion collisions is of femto-scopic in size, short-lived, and dynamically evolving with time, which limits the growth of the correlation length. Model calculations suggest that the finite size and finite time effects attained in high energy heavy-ion collisions limits the value of the ξ ≈ 2-3 fm [24]. This makes it extremely challenging to observe signals of criticality in experiments. Fortunately, model-based calculations suggest that the higher the order of the cumulants, the more is the sensitivity to the correlation length [25]. Further, it is expected that non-Gaussian features in net-particle multiplicity distributions will increase if the system freezes-out closer to the critical point. Hence, higher-order cumulants are favored as they could provide enhanced and observable critical signals. The study of cumulants also provides measurements that can be used to establish thermalization in heavy-ion collisions, to find the crossover temperature between normal nuclear matter and a deconfined phase called the quark gluon plasma, and set a scale for the phase diagram of QCD [26]. Cumulants up to the 4 th order (C n , n ≤ 4) of event-by-event distributions of net-charge, net-kaon and net-proton were measured by the STAR Collaboration in the phase I of Beam Energy Scan (BES) program at RHIC [27,28,29,30,31,32].
This review will focus on the most recent measurements on net-particle cumulants [27,28,31,32] and their implication in the hunt for the QCD critical point. Net-proton and net-kaon are used as proxies for net-baryons and net-strangeness, respectively, as charge-neutral particles are generally not identified in the experiments on an event-by-event basis. Theoretical calculations have shown that the net-proton number fluctuations could reliably reflect the baryon number susceptibility in the vicinity of critical point [33]. The prescription to link net-proton number fluctuations to that of net-baryon are discussed in Ref. [34].

Theoretical status of QCD critical point
First principle calculations for strongly coupled systems are, in general, very challenging. Theoretically, one would need to calculate the partition function of the QCD and obtain the singularities corresponding to the first-order phase transition line to locate the (T , µ B ) of the CP. This demands solving infinite-dimensional integrals, which are beyond our current capabilities, and hence one resorts to the numerical methods, i.e. the lattice-QCD (LQCD) simulations. At zero µ B , LQCD allows calculation of the equation of state of QCD as a function of T , indicating the nature of phase transition to be a smooth crossover. However, at finite µ B , which is the most suitable regime for the search of QCD critical point, LQCD suffers from the sign problem. The sign problem can be understood as follows: Typically, for any lattice computation one needs to evaluate the expectation value of an observable X, X(m v ) = , where M is the Dirac matrix in co-ordinate, colour, spin, flavour space for sea quarks of mass m s , S G is the gluonic action, and the observable X may contain fermion propagators of mass m v . The Det M for non zero µ is not positive definite, hence numerical methods of evaluation of the expectation values are difficult; this is commonly referred to as the sign problem. There are some suggested ways to overcome this issue. (i) Reweighting the partition function in the vicinity of transition temperature and µ = 0 [35], (ii) Taylor expansion of thermodynamic observables in µ/T about µ = 0 [36], and (iii) Choosing the chemical potential to be imaginary will make the fermionic determinant positive [37]. Due to these difficulties in lattice calculations in the finite or large µ B region, one requires QCD-based theories to make predictions on the presence and location of the critical point. Various theoretical calculations (LQCD and QCD-based models) on the location of critical point are shown in Fig. 3. The location of QCD critical point from theory calculations is scattered over region of µ B = 200-1100 MeV and T = 40-180 MeV [38,35,39,40,41,42,43,44,45,46,47,48,49,50,51,52]. While the accuracy of predictions for CP from the first principle lattice-QCD calculations worsens towards very large µ B , various model calculations wildly vary in their predictions. Therefore, experimental search of critical point is crucial to help establish it in the phase diagram.  [38,35,39,40,41,42,43,44,45,46,47,48,49,50,51,52]. Also shown on the top are the regions in µ B covered by various facilities and experimental programs.
One of the widely used theoretical calculations for the search of QCD critical point is with the linear sigma model [25]. A well-known behavior in critical phenomena is the divergence of response function at criticality due to sudden growth of the correlation length. At the CP, the exponents of correlation length on which various orders of thermodynamic response functions depend on are universal [53] and hence only rely on degrees of freedom in the theory and its symmetry. In view of the static critical phenomena, the QCD critical point falls into the Ising universality class [54]. The linear sigma model calculations qualitatively predict the universal critical behavior of the kurtosis of the order parameter fluctuations, i.e., of the σ-field, near the QCD critical point. Introducing coupling of the particles to the σ-field, allows the connection of fluctuations of the order parameter field, σ, to the fluctuations of the experimentally observable quantities, like the multiplicity distributions. The kurtosis of the order parameter, the coupling of protons to the σ-field, as well as the kurtosis of proton multiplicity distribution is discussed in Ref. [25]. To exploit the universality of critical phenomena, mapping of Ising model coordinates r and h (r is the reduced temperature and h is the external ordering field in Ising) to the T − µ B plane of QCD phase diagram is needed. This requires the knowledge of QCD equation of state with a critical point. Using information from lattice-QCD at µ B =0, and known features of thermodynamic variables near the QCD critical point, such a mapping is discussed in Ref. [25,55], where certain predictions about the behavior of QCD near the critical point was made. As shown in the Fig. 4, upon approaching the critical point from the crossover side of phase transition, the kurtosis of the order parameter is universally negative, and hence kurtosis of fluctuation observables like proton multiplicity distributions will be driven to take lower values than purely statistical fluctuations. Receding away from the critical point towards the region of first-order phase transition, the effect which drove the negative kurtosis of the order parameter previously will cause a higher value of kurtosis of proton multiplicity distributions as compared to baseline statistical fluctuations. Hence, non-monotonic energy dependence of the normalized kurtosis of proton multiplicity distributions with respect to baseline fluctuations would suggest the presence of a QCD critical point. The magnitude of the deviation from the baseline fluctuations will depend on the proximity of the freeze-out to the critical point. This model considers the critical fluctuations with the assumption of a static and infinite medium at equilibrium. However, one needs to keep in mind the finite size and time effect of the system created in heavy-ion collisions. Effects like initial geometry fluctuations, diffusion of conserved charge, the time evolution of fluctuations, jets, the role of conservation laws, resonances could also influence the fluctuation measurements in heavyion collisions. Nevertheless, the non-monotonic energy dependence of the kurtosis in the presence of a critical point has been confirmed by various other QCD-based model calculations like the Nambu-Jona-Lasinio (NJL) [56,57], Polyakov loop extended Quark Meson Model (PQM) [58,59,60]. Functional renormalisation group (FRG) approach to QCD also suggest non-monotonic collision energy dependence of fourth-order baryon number susceptibility(χ B 4 /χ B 2 ) due to sharpening of chiral crossover at large µ B , which may hint at a change in nature of phase transition at higher µ B [61].
Confirmation of a first-order phase transition would indeed validate the presence of a critical point. This is true, given that lattice-QCD calculations have confirmed the transition to be a crossover at µ B = 0 [8]. If the system formed in heavy-ion collisions is in the proximity of first-order phase transition, the multiplicity distributions could become bi-modal/two-component in nature, reflecting the contributions from two distinct phases. The two-component distribution has a very specific behavior of factorial cumulants, i.e.; they increase in magnitude with increasing order and alternate in sign [62,63,64]. This characteristic behavior is true even if one of the components is very small. This can be verified in the experimental measurements to search for first-order phase transition at finite µ B . Also, there has been no direct experimental evidence of crossover in the small µ B region so far. LQCD calculations have been extended to non-zero values of µ B using Taylor series expansion about vanishing µ B . Recent LQCD results for µ B < 110 MeV, suggests negative sign of ratio of fifth-to-first and sixth-to-second order baryon susceptibility [65]. The FRG model calculations also give the negative sign for the two ratios across a large range of µ B which is accessed in gold nuclei collisions at √ s N N = 7.7 -200 GeV  [25]. The y axis is the reduced temperature, or the difference of temperature in the Ising model from the critical temperature and the x axis reflects the external ordering field which results in breaking of symmetry in the Ising model.
in the STAR experiment at RHIC [61]. A PQM model study also suggests that for freeze-out close to the chiral phase transition temperature, the ratio of sixth-to-second order baryon susceptibility remains negative [58]. The nature of phase transition in these calculations is a crossover. Hence, higher-order cumulant measurements can serve as a litmus test to identify the nature of phase transition.

Experimental search for QCD critical point
This section deals with various aspects related to the experimental search for the QCD critical point. First, the observables employed in the search for the QCD critical point are presented. Baselines from various non-critical physics, like the thermal statistical fluctuations, baryon number conservation, and resonances, are also discussed. There are several complexities involved in experimental measurements of the number fluctuations like the volume fluctuations, efficiency correction for the high order cumulants, and centrality resolution. The origin of these experimental challenges and how they are addressed while performing the measurements are dealt with in the subsequent subsections. Finally, the results from the experimental measurements are presented and compared with various model calculations to understand their physics implications.

Observables
Cumulants quantify the traits of a distribution. The cumulants of a distribution up to the sixth order are defined as follows: where N is the observable whose distribution is being considered, which for our case is the net-particle number in a single event and N is the average of the N from all the events. The difference of N from N is given by δN = N − N . The first and second order cumulants are the well known mean (M) and variance (σ 2 ) of a distribution whereas the third and fourth order cumulants reflect the skewness (S) and kurtosis (κ) of a distribution, respectively. Cumulants can be related with conserved charge number susceptibilities for a system in thermal equilibrium. In LQCD and hadron resonance gas (HRG) model calculations, using the grand canonical ensemble construct, the n th order number susceptibility (χ q n ) associated to conserved charges, q (q = B, Q, S) can be calculated by taking derivative of dimensionless pressure (P/T 4 ) with respect to reduced chemical potential (µ q /T ) corresponding to the conserved charge q and is given as, [66,26] Here P and T represents the system pressure and temperature, respectively. Pressure can be written in terms of logarithm of partition function (Z) as follows: One can relate the cumulants of the conserved charge multiplicity distribution to the corresponding number susceptibilities as follows.
As can be seen from the above equation, all the cumulant orders directly depend on system volume (V ). Constructing the ratio of cumulants eliminates the trivial system volume dependence and thus allows for a direct comparison of the experimental measurements to the thermodynamic susceptibilities calculated in lattice-QCD, HRG, and various QCD-based models, although subject to some caveats. The experimental measurements are carried out within specific kinematic phase space, and probe the conserved charge quantities via forming proxies from the detected charged particles. The theoretical calculations like those from the lattice-QCD and QCD-based models instead deal with conserved charge quantities and do not usually involve any phase space cuts. These limitations should be kept in mind when comparing the experimentally measured cumulant ratios to the susceptibility ratios obtained from theoretical calculations presented throughout this review article. In the existing literature, both the ratios of cumulants and the product of moments of the distribution are used. The most widely used products of moments in the literature are the σ 2 /M, Sσ, and κσ 2 . The exact relation of these moment products to the cumulants is listed below.
From cumulants, one can construct the factorial cumulants, also known as the integrated correlation functions. The factorial cumulants (κ n ) are related to cumulants (C n ) as follows: In addition, extensive studies on various other fluctuation observables have been done in the existing literature. Dynamical K/π, p/π and K/p fluctuations in terms of an observable which is robust to detector efficiency effects called the ν dyn have been reported by experiments [67]. Another such robust observable is the balance function, which measures probability of observing net-particle separated in pseudo-rapidity. Narrowing of balance function was observed by the experiments [68], which is consistent with the concept of delayed hadronization of a deconfined QGP. Measurement of p T dynamical correlations have also been performed by experiments [69]. Monotonic scaling of the measurements on dynamical correlations has been observed so far.

Non-critical point contributions and baselines
This section discusses various non-critical baseline calculations. Although a non-monotonic collision energy dependence has been suggested as being the signature of the presence of QCD critical point, there could be several other non-critical effects causing appreciable correlations in the experimental measurements. For example, the effect of unavoidable physical effects like the purely statistical fluctuations, baryon stopping at lower energies, baryon number conservation, and resonances [70,71,72,73,118], needs to be studied. In this regard, in addition to statistical baselines, two widely studied models, which do not incorporate any physics of QCD critical point or QCD phase transition, are discussed, namely the HRG [75,76] and the UrQMD [77].

Statistical baseline
If the fluctuations of net-particle multiplicity distributions are purely Gaussian in nature, the third and all higher-order cumulants become zero. Hence, the presence of non-zero third, fourth, and all higher-order cumulant is the first sign to take note of while probing critical behavior of cumulant. One of the widely used statistical baselines is the Skellam baseline. A Skellam distribution is the resulting distribution of difference of two variates that follow independent Poisson distributions. For example, if the protons and anti-protons follow independent Poisson distributions, the distribution of net-protons is Skellam. The Poisson distribution has the unique property of having the same value for all the order of cumulants. Cumulants, in general, follow additive property, i.e., given two independent random variables X and Y , cumulants of a new variate Z, constructed by addition of X and Y (Z = X + (−1) m Y , where m can be any integer), could be obtained by simply adding the cumulants of the input variates (C Z n = C X n + (−1) n C Y n ). Using the previous two statements, it can be easily shown that the ratio of odd-to-odd and even-to-even cumulants calculated from Skellam expectations are unity. Similarly, one can also assume binomial and negative binomial distribution for the constituent variates of the net-particle distributions for statistical baseline. The table 1 tabulates the values of statistical baseline for cumulants of net-particle distributions (Z) up to the sixth-order, when the constituent variates of the distributions (X and Y : Z = X -Y ) follow independent Poisson, Gaussian, Binomial, and Negative Binomial distributions(NBD). In the table, the baseline expectations for Binomial and NBD distributions are expressed in terms of their cumulants. The cumulants of a Binomial distribution up to sixth-order can be obtained in terms of its parameters; the number of trials (n), and probability of success (p) as follows: Similarly, cumulants of an NBD distributions up to the sixth order can be expressed in terms of its parameter; the number of failures (r), and probability of success (p) as follows: , a Gaussian has two parameters: the mean (µ) and width (σ), Binomial has two parameters: number of trials (n) and probability of success (p), and NBD also has two parameters: number of failures (r) and probability of success (p), respectively. Cumulant Poisson Gaussian Binomial NBD With the knowledge of experimental measurement of first and second-order cumulant of the constituent distribution, one can extract the parameters for the underlying statistical distributions considered above for statistical baseline and make predictions for the higher-order cumulants of the net-particle distributions. For Poisson distribution, only the experimental measurement of mean is required as input as it has only one parameter.

Hadron resonance gas model
In the ideal or non-interacting hadron resonance gas (HRG) model, the degrees of freedom is comprised of point-like hadrons and resonances which are in thermal equilibrium. The logarithm of the partition function of a hadron resonance gas in the grand canonical ensemble (GCE) is given as where the summation is over all hadrons and resonances. Further, for each particle, where, T is the temperature, V is the system volume, E i is the energy, g i is the degeneracy factor and µ i = B i µ B + S i µ S + Q i µ Q is the chemical potential of the i th particle. Here in the last expression, B i , S i and Q i are the baryon, electric charge, and strangeness number of the i th particle, with corresponding chemical potentials µ B , µ Q , and µ S , respectively. The ± sign are for baryons and mesons, respectively. Once the partition function of the hadron resonance gas is known, using the eqns 7 and 8, the thermodynamic susceptibilities can be obtained. At higher temperatures and large chemical potential values, the ideal gas assumption may no longer be valid. Van-der-Waals(VDW)-type interactions can be incorporated into the HRG to reflect the qualitative feature of strongly interacting gas of hadrons. In the excluded-volume (EV) approach to HRG, Van-der-Waal type repulsive interactions are introduced by relaxing the argument of the point-like size of hadrons and considering the geometrical size of hadrons. The Van-der-Waals equation of state is given as follows: Here, P , V , and N are the pressure, volume and number of particles of the system, respectively. The parameters a and b are the Van-der-Waals parameters which describe the attractive and repulsive interactions, respectively. One can also include both attractive and repulsive Van-der-Waals interactions to the HRG model framework. Susceptibility calculations using varieties of the HRG models in the existing literature are discussed in Ref. [75,76,78,79,80,81,82]. The most commonly used variants of HRG (both ideal and interacting) employs the Grand Canonical Ensemble (GCE), which results in the conservation of quantum number on an average. However, in the low center-of-mass energies for heavy-ion collisions, the canonical ensemble (CE) approach becomes more suitable, which conserves the exact quantum number as opposed to the GCE construct [83]. The fluctuation in net-proton is constrained to the net-baryon number in the full phase space, which itself is a conserved quantity. Using the information on the mean number of accepted protons and antiprotons within the kinematic acceptance, fluctuation of net-proton in the selected kinematic acceptance is modelled [83].

UrQMD model
The Ultra-relativistic Quantum Molecular Dynamics (UrQMD) model is a microscopic transport model [77]. The model includes transportation of various degrees of freedom (e.g. baryons and mesons) and the production of new particles and their interaction. The production of particles occurs via fragmentation of strings made of valence quarks of the original colliding hadrons, resonance excitation, and decays. There is no quark-hadron phase transition implemented in this model. To study the effect of baryon number conservation, baryon stopping at lower energies, and resonances on the fluctuation measurements, UrQMD calculations serve as a suitable baseline for the collision energies at RHIC.

Experimental challenges
Although obtaining the number of particles in an event may appear simple, several experimental techniques are employed in the measurements of net-particle fluctuations. Event statistics required for proper estimation of cumulants, purity of the charged particles selected, centrality/impact parameter estimates, removing self-correlation between particles used to define centrality and net-particle fluctuation, volume fluctuations and ways to suppress the effect, correction of cumulants for finite detector efficiency, estimation of uncertainties, are some of the topics which are of prime importance in fluctuation measurements. Current subsection discusses these aspects of the measurements.

Event statistics and statistical uncertainties
As experiments run for a finite amount of time and collect a limited amount of data, estimates should be made a priori on how much event statistics are sufficient enough to measure an observable of interest. The Ref. [84] performs such a model study and presents the minimum event statistics needed for estimation of cumulants (C n ; n < 8) in order to detect a signal (related to phase transition effects or critical point) of 5% and 10% above the statistical baseline. The event-by-event net-proton distribution constructed as the difference of the number of protons and anti-protons was simulated using two types of underlying distributions, Poisson and Binomial. The input parameters are taken as the mean and width of (anti-)proton distributions in 0-5% central Au+Au collisions at √ s N N = 62.4 and 200 GeV measured by the STAR detector [30]. Minimum statistics required to estimate the cumulants of netproton distribution with a precision of 5% with Poisson being the nature of underlying distributions for protons and anti-protons are shown in Fig. 5. The limits are shown for two sets of input parameters determined from the STAR experiment. Among various orders, we discuss here the statistics needed for C 4 and C 6 measurements. Estimation of C 4 (C 6 ) within a precision of 5% requires 1.2 million (1.5 billion) event statistics with input parameters taken from √ s N N = 200 GeV while they are found to be 1.6 million (1.5 billion) with input parameters from √ s N N = 62.4 GeV. Similar study done with binomial assumption for protons and anti-protons distribution in the simulation fixing the inputs from √ s N N = 200 GeV suggests measurement of C 4 (C 6 ) at a level of 5% precision demands 1.3 million (1.8 billion) events [84]. Taking the inputs from √ s N N = 62.4 GeV also gives similar event statistic limits.
One of the important aspects of any measurement is the proper estimation of uncertainties or errors on them. We discuss in this subsection how the statistical uncertainties on the fluctuation measurements are calculated. The limited event statistics available for measurements demands a careful estimation of statistical uncertainties. Statistical uncertainties on the cumulants are usually obtained using the Delta theorem method and Bootstrap method. While the Delta theorem method is an analytical method of standard error propagation, the Bootstrap method is a computer-intensive re-sampling method.
Using the Delta theorem, statistical uncertainties on cumulants and moments have been computed [85,86] and have been quite extensively used in the analysis of higher-order cumulants in heavy-ion collision experiments. This error estimation method has at its foundation an important theorem on the convergence of probability known as the central limit theorem (CLT), which goes as follows. Central Limit Theorem: Suppose X 1 , X 2 , .., X n is a collection of random variables that are independent and identically distributed (i.i.d.) with E[X i ] = µ and Var[X i ] = σ 2 , then in the limit n approaches infinity, the random variable √ n(S n − µ) converge to a normal distribution N(0, σ 2 ), where S n = (X 1 + X 2 + .. + X n )/n. In other words, for large value of n, S n approximately follows a normal distribution N(µ, σ 2 /n). Simply put, the CLT theorem states that several important random variables and estimators are asymptotically normal. The Delta theorem method allows approximation of the asymptotic behavior of functions over a random variable if the random variable itself is asymptotically normal [87]. The mathematical statement of delta theorem is given as: (Delta theorem) Let T n be a sequence of statistics such that Let g be a real function which at least is differentiable at θ with g ′ (θ) = 0. Then, Here, the d − → sign represents convergence over the distribution. It is worth mentioning here that the estimation of error using the Delta theorem involves the parameters of the population (e.g. σ and µ in the CLT definition), which is inaccessible, as one only has a fixed sample to work with. Hence, the parameters of the population are estimated using the sample itself. Delta theorem can be easily extended to multivariate case. Then, where, X = {X n } is sequence of random vectors which are normally distributed; X d − → N(θ, Σ/n), θ is a constant vector (mean vector), Σ is the covariance matrix, and D is the Jacobian of g. For estimating error on the cumulants, we will make use of one more theorem employed in convergence of a sequence of sample moments, which actually comes as a consequence of multivariate central limit theorem.
where, Σ is a (k − 1) × (k − 1) covariance matrix, with elements, Using the above theorem and putting i = j = m (m >= 2) in the eqn. 38, the variances of m th order central moments can be obtained. Using the CLT and the above theorem, error on the cumulants up to the third order are obtain as follows: To avoid the use of square root symbol repetitively, expressions for variances of the cumulants are presented in eqn. 39. The statistical uncertainties/errors on the cumulants will simply be the square roots of the variances. Obtaining error on the fourth order cumulant will require the use of Delta theorem method given in eqn. 36. Taking g(µ 2 , µ 4 ) = C 4 = µ 4 − 3µ 2 2 , Since the central moments µ 2 and µ 4 are involved in the transformation, the elements of the covariance matrix are Σ 11 = V ar(µ 2 ), Σ 12 = Σ 21 = Cov(µ 2 , µ 4 ), and Σ 22 = V ar(µ 4 ). The variances and covariances of the central moments can be calculated using eqn. 38. Then, Using the above equation, the variance on C 4 is as follows: With increase in the order of cumulant, the formula for uncertainties on the cumulant becomes more tedious. Involving the detection efficiency while calculating the statistical uncertainties complexifies the formulae further [86]. It is seen that larger the width of a distribution and/or smaller the detection efficiency, the larger is the statistical uncertainties on the cumulants. As an alternative to the complicated error formulae, especially for the higher-order cumulants, it is quite helpful sometimes to use re-sampling methods over the analytical delta theorem method for estimation of statistical uncertainties on cumulants.
(Bootstrap Method) Bootstrap method is a re-sampling method to estimate statistical error and thus does not involve the complexities of the standard error propagation method. It uses Monte-Carlo algorithm to estimate the statistical error on a statistic by creating many bootstrap samples, where each bootstrap sample comprises of elements randomly drawn from the original sample with replacement. The statistic is calculated for each of these bootstrap samples and the sampling variance of the statistic from these bootstrap samples is the error on it [88,89,84]. To start with, let X be a random sample with n data points representing the experimental data set randomly drawn from an unknown parent distribution. Let the estimator be denoted byê, which could be any statistic such as mean or variance, whose standard error we intend to find. The prescription to estimate standard error using the bootstrap method is as follows: 1. Construct B number of independent bootstrap samples X * 1 , X * 2 , X * 3 , ..., X * B , each consisting of n data points randomly drawn from the random sample X with replacement. 2. Evaluate the estimator in each of these bootstrap samples, 3. The sampling variance of the estimator is given as follows.
. The sufficient enough value of B for an accurate estimation of error within the bootstrap method varies from case to case depending upon the initial sample size. However, in general, the larger value of B estimates the error better. Using a Monte-Carlo procedure, both the methods were subjected to verification of Gaussian nature of the statistical error and found to satisfy the criteria [84]. Comparison of error estimated by both the methods on fourth order net-proton cumulant C 4 for Au+Au collisions at √ s N N = 54.4 GeV from the UrQMD model as a function of N part are shown in Fig. 6. The errors estimated from both the methods show good agreement with each other.

Particle identification
The results presented in this review dominantly cover the measurements from the STAR experiment, hence we discuss the particle identification effects using it as a typical example. The STAR detector at RHIC using two sub-detectors Time Projection Chamber (TPC) and Time of Flight (TOF) ensures the high purity of the identified particles [90]. TPC measures the momentum, trajectory and ionization energy loss per unit length (dE/dx) of charged particles. Only the particles with transverse momentum greater than 150 MeV/c are recorded. The dE/dx measurement is used for particle identification (PID) by TPC. Figure 7 shows average dE/dx as a function of p/q (momentum/charge) in Au+Au collision at √ s N N = 39 GeV at RHIC-STAR [91]. A variable called the nσ X is used for selection of identified particles, which can be constructed from the dE/dx information, as shown below, Figure 7: (dE/dx) of charged tracks measured by TPC as a function of p/q (momentum/charge) in Au+Au collision at √ s N N = 39 GeV in STAR [91]. The various bands represents different charged particles, namely the protons (p), anti-protons (p), kaons (K + , K − ) and pions (π + , π − ). The curved lines in black are the Bichsel expectation for various charged particles [92].
where, (dE/dx) measured is the ionisation energy loss of a charged particle measured by TPC and (dE/dx) theory is the expectation from Bichsel formula [92] for the charged particle. σ X is the dE/dx resolution of the TPC. The particle identification with |nσ X | < 2 is found to give a good level of purity (> 90% for kaons and > 98% for protons [28,30]) of the identified charged particle in the momentum range where only TPC is used for fluctuation measurements. For the lower momentum region, the bands corresponding to various charged particles are well separated while at the larger momentum, the bands can be seen to overlap with each other and thus the purity deteriorates. To identify charged particles in the higher momentum region, the TOF detector is used. The TOF detector measures flight time (t) of charged tracks from primary vertex (interaction point where collision between two heavy-ion happens) to the detector. The square of mass of charged particle which requires flight time information from TOF along with information on momentum (p) and charge (q) from TPC, is used for particle identification. The mass squared (m 2 ) as a function of p × q for Au+Au collisions recorded by STAR detector at √ s N N = 54.4 GeV is shown in Fig. 8. If L is the path length, v is the velocity of the charged particles, then β can be written as β = v/c = L/ct. The m 2 can then be written as, The above equation suggests for a given momentum (p), charged particles of larger mass (m) will lead to larger flight time (t) and vice-versa. The m 2 measurement by TOF allows for PID with high degree of purity (> 90% for kaons and > 97% for protons [28,31]) in the higher momentum region.

Linking measurements to conserved quantities
As discussed in the previous section, only charged particles could be measured using the TPC detector.
Measuring net-baryon number fluctuations poses a challenge as neutrons, which are also produced in large numbers in experiments and carry a unit baryon number, are not identified. Interestingly, it was found that net-proton number fluctuations could reliably reflect the baryon number susceptibility near the critical point [33]. Theoretical works to link net-proton cumulants to those of net-baryon have been presented in Ref. [93,34]. This conversion of cumulants from those of net-proton to net-baryon relies on the fact that isospin of nucleons in the final state are randomized and uncorrelated [93]. Correlations between the isospins of nucleons in the final state are found to be negligible over a wide range of collision energy. In central heavy-ion collisions at large energies, the isospin density becomes negligibly small as most of the initial isospin density is absorbed by pions which are abundantly produced. Then the conversion of cumulants of net-proton (C net np ) to those of net-baryon (C net nB ) gets simplified and are given as follows [93,34]: where C tot np is the n th order cumulant of proton+anti-proton distribution and C net,tot np,mp are the mixed cumulants of net-proton and proton+anti-proton distribution of (n + m) th order.

Centrality determination
In heavy-ion collisions, the impact parameter (b) of collisions is defined by the distance between the geometrical centers of the colliding nuclei in the plane transverse to their direction. A schematic diagram of the geometry of heavy-ion collisions is shown in Fig. 9. As the impact parameter is not measurable in experiments, the charged particle multiplicity of events can be used to define collision centrality as they reflect the initial geometry of the colliding nuclei. To obtain the geometrical quantities like the impact parameter, the number of participant nucleons (N part ) and the number of binary collisions (N coll ), one has to rely on model calculations which involves mapping the charged particle multiplicity from experiment to that of a simulated one [94,95,96]. The widely used model for determining these quantities is the Glauber Monte Carlo (Glauber MC) model [96]. In Glauber model, nucleons in a nuclei are distributed following the Wood-Saxon density distribution. The nuclei are then assigned a random impact parameter and translated to collide. The distance between the centers of nucleons of each nucleus is calculated. If this distance is found to be less than σ N N inel /π, where σ N N inel is the inelastic nucleonnucleon cross section, the nucleons have undergone a binary collision. Nucleons that have undergone at least one binary collision are said to have participated in the reaction. This procedure is followed to obtain N part , N coll event-by-event. A class of events belonging to a given centrality (say 0-5%) is obtained by the collection of events representing the fraction of the total cross-section. Glauber model can only provide the information on N part , N coll and impact parameter b, but for mapping with the experimental data, one has to use a combination of the Glauber model and a particle production model called the Two-Component Model [97]. Assuming that particle production is governed by contribution from hard component (proportional to N coll ) and soft component (proportional to N part ), it simulates the multiplicity distribution according to the eqn 51. 0%-5% 5%-10% 10%-20% 20%-30% 30%-50% Figure 10: An illustration of the total inclusive charged particles multiplicity N ch and number of participant (N part ) and impact parameter (b) in heavy-ion collisions [96]. Figure taken from Ref. [96].
where the quantity n pp is the multiplicity from a p + p collision of the same center-of-mass energy.
x is the contribution of hard component. n pp is drawn from a negative binomial distribution with parameters NBD (n pp ; n pp ,k), where n pp is the average multiplicity in the p + p collision and k is the width of the distribution. The simulated charged multiplicity distribution is then fitted with the corresponding multiplicity distribution from the experiment. n pp , k, and x are kept as free parameters in the fitting to experimentally measured charged particle multiplicity distribution. Initial guess values for x could be set close to those determined in different experiments [98,99] for faster convergence of fitting. Figure 10 illustrates collision centrality definition in heavy-ion collisions by comparing the charge particle multiplicities with Glauber MC + Two-component model simulation. Centrality classes are classified for the simulated distribution. N part , N coll in those centrality classes are calculated.
In addition to using the number of charged particles to define centrality, one can also use detectors different from the ones used for particle identification, such as the Zero Degree Calorimeters (ZDC) and the Event Plane Detector (EPD). These detectors measure the energy deposition in the forward rapidity region and thus could help eliminate any possible correlations between the cumulants measured at mid-rapidity and the centrality.

Volume fluctuations
The impact parameter of the collisions could change event-by-event. The fluctuation of the impact parameter and thus the participant nucleons gives rise to fluctuations in the initial system volume created. This effect is called the volume fluctuations (VF), which could result in additional nondynamical fluctuations in the experimental measurements. Studies on the effect of VF on net-particle cumulants suggest that there is an artificial enhancement of the cumulants due to VF. The methods to correct the cumulants for this effect within the limitation of experimentally defined collision centrality will be discussed in this sub-section.
Centrality Bin Width Correction (CBWC): Any given centrality, say 0 − 5%, corresponds to a range of impact parameters. Selection of narrow centrality bins helps to get rid of inherent fluctuations within each centrality bin. For the charged particle multiplicity distribution, the smallest bin-width can be of unity-width. The correction to cumulants of net-particle distributions associated with accounting for the finite centrality bin width is called the centrality bin width (CBW) correction [100]. For a given centrality that spans over a range of charged particle multiplicity, cumulants in each multiplicity bin are evaluated. A weighted average of cumulants is made with the events in each multiplicity bin as weights to quote the CBW corrected value of the cumulant for that centrality. The CBW corrected n th order cumulant in a certain centrality is given as where r , sums over all multiplicity bins in that centrality. n r is the number of events in the r th multiplicity bin. N is the total number of events in the centrality considered. The error on the CBW corrected cumulants is evaluated using standard error propagation method to be: where σ Cn,r is the statistical error on n th cumulant in r th multiplicity bin.
Volume Fluctuation Correction (VFC): Volume fluctuation correction [71,101] is another method to suppress volume fluctuations effects on cumulant measurement. The correction procedure assumes each N part value to be an independent source of particle production. Net-particle numbers (∆N) for a given centrality class are constructed as sum of net-particles (∆n) from each N part source. As the N part itself fluctuates in a given centrality, the cumulants of net-particle distributions get contribution due to the N part fluctuation. The VFC method requires the estimation of cumulants of N part fluctuation in addition to those of net-particle fluctuations. As the N part is inaccessible in experiments, particle production models which have information on geometrical quantities are used in VFC. The correction prescribed by VFC method for cumulants up to the 4 th order is as follows [71,101]: The above equations suggests that starting from second order cumulants of net-particle distributions, the N part fluctuation encoded in their higher order cumulants (C n (N p ), n ≥ 2) also contributes. In the limit ∆n − → 0, which is realized at very large collision energies, the N part fluctuation contributions up to the third order net-particle cumulants vanish, but there could be sizable contributions for the fourth order cumulant.
Comparison of these two methods for the net-proton C 4 /C 2 at three RHIC energies, 7.7, 19.6 and Figure 11: 19.6, 62.4 GeV [32]. Results from CBWC (red circles) and VFC (blue squares) method to correct for volume fluctuations are presented.
62.4 GeV for 0-5% centrality are shown in Fig 11. The correction factors for the VFC method was obtained from Glauber+Two component model [96]. The results of CBWC and VFC are found to be consistent with each other within uncertainties. While the CBWC is data-driven correction method, the VFC is model-dependent. While performing CBWC, one needs to note that even in the smallest centrality bin limit of unit multiplicity, participant fluctuation could still persist, although arguably small. The VFC, on the other hand, requires the assumption of independent particle production from sources. Model studies by authors of [102] suggest that the VFC method does not work well with the UrQMD model, which is attributed to the breakdown of the assumption of independent particle production in the UrQMD.

Centrality resolution and self-correlation effects
Experimentally, while defining centrality using charged-particle multiplicity, one needs to carefully study two important aspects: centrality resolution and self-correlation effect. For the CBWC method to sufficiently suppress the volume fluctuations, the resolution of the centrality estimator should be high.
Increasing the charged particle number in centrality definition (for example, by increasing the pseudorapidity (η) acceptance) improves the centrality resolution. Figure 12 presents a UrQMD model study on the dependence of net-proton Sσ and κσ 2 on different centrality definitions obtained using varying η acceptance [100]. With a smaller η acceptance used for centrality definition, an artificial rise is found for the two moment products at all the energies, which is suppressed when larger η acceptance are used to define centrality. However, that could lead to self-correlation effects if the charged particle acceptance used for centrality definition overlaps with the acceptance taken for net-particle selection for fluctuation measurements. Self-correlation effect arises when the charged particles which are selected for fluctuation  Self-correlation effect could also arise due to decay of resonances where one or more of the decay daughters falls into analysis acceptance, and the other daughters are taken in centrality definition. Weak decays of lambda (Λ), cascade (Ξ) and omega (Ω) baryons could decay into protons (anti-protons) and charged pions or kaons. Self-correlation could arise if protons from such decays are selected for analysis and pions and kaons are included in centrality selection. There are also other resonances ∆ ++ and ∆ 0 which could contribute to the effect. Although keeping a tight cut on a track's distance of closest approach (DCA) to the primary vertex suppresses the background from decays and resonances in particle identification to a good extent, but it cannot be fully eliminated. For example, the MC simulation studies done by the STAR experiment with HIJING [103] as input indicate that about 81% of the lambda-decay protons were reconstructed with DCA < 3 cm and 65% with DCA < 1 cm [104]. As the fluctuation measurements at STAR consider a criterion of DCA < 1 cm for particle identification, the study of self-correlation arising out of such decays becomes important. A HIJING model based study to understand ∆ and weak decays of strange baryon induced self-correlation effects on net-proton cumulant ratio C 4 /C 2 is shown in Fig. 13. The C 4 /C 2 at √ s N N = 7.7 and 200 GeV from various centrality definitions aimed at removing such correlations are found to be consistent with each other, suggesting that the self-correlation effects from resonance decay are negligible. For the discussed fluctuation measurement results in this review, the effect of pseudo-rapidity range for centrality determination has been understood and optimized using model calculations in order to maximize resolution and suppress self-correlations [31].

Pileup rejection
The TPC detectors typically have a drift time in the orders of tens of µs. If this time duration is longer than the bunch crossing rate of the colliding ions in heavy-ion experiments, then there are possibilities that all the charged tracks from an event are not collected, and the subsequent events are triggered. In such scenarios, where successive events are reconstructed as a single event, there is a piling up of charged particles, and the event is called a pileup event. If such events are sufficiently large in number, they can alter the event-by-event net-particle distributions and consequently their cumulants. Detectors that have a faster detection time are used to remove pileup events, for example, the TOF detector at STAR experiment and the START detector at HADES experiment. The correlation between the charged particle multiplicity registered at the TOF and the TPC is examined to remove the background coming from pileup effect. It utilizes the fact that the TOF is a much faster detector than the TPC with a timing resolution of 100ps and thus sees less track fraction from previous events. The high interaction rates in the fixed target (FXT) configuration of experiments result in large pileup effects compared to experiments run in collider mode. Removing pileup events then becomes a challenging task. The excellent timing properties of the START detector at HADES help reduce the contamination of identified proton yields by pileup to a probability of ≤ 3 × 10 −5 in the proton number fluctuation studies in Au+Au collisions collected in the FXT mode at √ s N N = 2.4 GeV [105]. The recent measurements of proton cumulants in Au+Au collisions at √ s N N = 3 GeV (FXT) by the STAR experiment [106], applied the correction of cumulants for pileup effect suggested in Ref. [107,108]. The correction deals with two distributions: one is the multiplicity distribution of produced particle (m) that is considered for centrality definition, and the other is the distribution of the particle of interest (N) whose fluctuations are to be calculated. The observed particle number distributions P m (N) at some value of multiplicity m is believed to have a contribution from true or single collision events P t m (N) and pileup events P pu m (N). Mathematically, P m (N) can be written as: where α m is the contribution of pileup events to the observed particle number at multiplicity m. The correction procedure assumes pileup events as the superposition of two independent single-collision events as suggested by the eqn. 60 [107].
Here, P sub i,j is the probability distribution of N in a pileup event which is constituted by two single collision events having multiplicity i and j. The combined multiplicity of the two events constituting the pileup, i.e., i + j must be equal to m. w i,j , here, refers to the probability of observing a pileup event at multiplicity m, which is constituted by two single events with multiplicity i and j.
The true cumulants in each multiplicity bin (cumulants of P t m (N)), then, can be obtained recursively from lower multiplicity bins applying the correction formulae prescribed in Ref. [107]. There are two parameters involved in the correction: α m and w i,j , which need to be determined from simulations. The correction's accuracy depends on how well the simulation mimics the real experimental conditions. Studies by the authors of Ref. [108] suggest that the parameters involved in the correction formula can be obtained by a model-independent unfolding approach. The whole formalism of the correction can be extended to situations where pileup can also be caused due to more than two single collision events, although the probability of such pileup events will be very small.

Efficiency correction
Binomial efficiency correction: The detection of particles in a detector can be modeled assuming Binomial sampling B ǫ,N (n) where ǫ (probability of success) is the efficiency parameter, and N (number of trials) is the total number of charged particles incident on the detector. The assumption of Binomial sampling has been tested in several works from the STAR experiment and is shown to be a good working assumption for event-by-event detection of charged particles [28,32,109]. The true multiplicity distribution (P (N)) gets modified due to detector efficiency. The measured multiplicity distribution is a convolution of true multiplicity distribution and distribution for detector response, which is assumed to be binomial. Individual particles are detected with a probability ǫ, which is assumed to be independent for different particles. Denoting the number of detected particles as n and the distribution of measured multiplicity asP (n), the relation between the measured and true distribution could be given as, As the experiments are interested in fluctuations of net-particle (n = n net = n 1 − n 2 ), the bivariate case of the above relation are used to obtain formulae for efficiency correction. Then, using the factorial cumulants generating function, the cumulants of the true net-particle distribution can be obtained using the cumulants of the measured distributions as follows [110]: here, n o =< n 1 + n 2 >. C n and c n are the cumulants of true and measured net-particle multiplicity distributions, respectively. f i,j are the factorial cumulants of the measured multiplicity distributions. This modeling can be extended to multivariate case where there could be different detection efficiency for different particle species and momentum ranges as is the case in most of the experimental measurements. The correction of cumulants for detector efficiency then gets more involved. A detailed discussion on the efficiency correction for many efficiency bins can be found in Ref. [86,111,112]. The complications and computing time involved in efficiency correction with several efficiency bins reduces as compared to the approach in Ref. [86], if cumulants and mixed cumulants are obtained for an observable q (r,s) defined as, where, n i is the number of particles in i th efficiency bin and a i , ǫ i represents the quantum charge number, efficiency of particles in i th bin, respectively. M is the total number of efficiency bins. The efficiency corrected cumulants in terms of the cumulants and mixed cumulant of the above observable are presented in Ref. [111]. Here, in the eqn. 66, ǫ i is the average efficiency of particle in a bin which requires information of the particle spectrum for that bin. A simplified procedure called the track-bytrack efficiency correction eliminates this aspect of correction [112]. If M → ∞ in eqn. 66, i.e., there are infinite number of efficiency bins and q (r,s) can be expressed as follows.
With M → ∞, each efficiency bin could now contain up to one particle. For n i = 0, there is no contribution to q (r,s) and one can rewrite eqn. 67 as, where n tot is the total number of particles in one event. As the efficiency correction using eqn. 68 only demands track-by-track efficiency, analytical formula of efficiency with respect to a track variable like p T could be used. The particle yield information within each efficiency bin is no longer required.
Unfolding approach: The analytical formula for efficiency correction hinges on the fact that the detector response is binomial. In very high-energy head on collisions, a large number of charged particles are produced. Such dense charged particle environment may result in non-binomial effects in detector response. In the presence of non-binomial effects, the detection probability may not be independent for each particle which could result in unwanted correlation in the measurements. The correction of cumulants in such situations could be done by a method called unfolding [113,114,115]. Let the true particle distribution P (N) be linked to the measured particle distributionP (n) by the equation,P which is identical to eqn. 61 except the detector response given by T (N, n) is no longer assumed to be binomial. What unfolding does is to invert the detector response matrix and apply it to the distribution at the measured level to obtain the true net-particle distribution. The equation below depicts the procedure.
Unfolding method, in general, demands rigorous computation and resources. It has been applied by the STAR collaboration as a cross-check to get back cumulants of true distributions for central Au+Au collisions at √ s N N = 200 GeV [31]. It was found to give values consistent with the cumulants from the binomial efficiency correction procedure within uncertainties.

Systematic uncertainties
Systematic uncertainties are broadly classified under four categories (a) those arising due to the choice of criteria for selection of charged particles, (b) due to specific particle identification method, (c) those associated with the estimation of efficiency of detectors, and (d) due to use of specific event selection criteria. Typically they are estimated by varying the following requirements for charged particle tracks: DCA, track quality reflected by the number of fit points used in track reconstruction algorithm in the TPC, number of fit points associated with the dE/dx measurements, dE/dx and m 2 for charged particle identification [27,28,30,32]. A ± 5% systematic uncertainty associated with determining the efficiency is also considered [91]. All of the different sources of systematic uncertainty are added in quadrature to obtain the final systematic uncertainties on the cumulants and their ratios. The typical value of statistical and systematic uncertainties on net-proton cumulant ratios C 4 /C 2 in most central 0-5% Au+Au collisions at √ s N N = 39 GeV is 20% and 18%, respectively [30].

Experimental results
In this section, the most recent results on the net-particle cumulant measurements from the Au+Au collisions recorded by the STAR detector [90]  GeV were analysed as part of the phase-I of the Beam Energy Scan (BES) program at RHIC [91,116]. Table 2 shows the total event statistics for the nine collision energies along with corresponding T and µ B at chemical freeze-out. Time-Projection-Chamber (TPC) and Time-of-Flight (TOF) detectors are used to select charged particles within p T range 0.2 -2.0 GeV/c. The measurement are carried out at mid-rapidity (|y| <0.5). The collision centrality is determined from the charged particle multiplicity excluding the particles of interest to avoid the self-correlation effect [117]. In order to suppress the volume fluctuation effects, centrality bin width correction is applied to the measurement of the cumulants [100]. Cumulants are corrected for acceptance and detection inefficiencies, with the assumption that the distribution of the detector response is binomial [86,111]. For estimation of statistical uncertainties of cumulants and their ratios, both Delta theorem method and Bootstrap are used [85,84]. Systematic uncertainties of the C n 's are estimated varying the tracking efficiency, track selection, and particle identification criteria. In the following subsections, the collision energy dependence of cumulants and cumulant ratios and comparison to various theoretical and model calculations are discussed.

Uncorrected net-proton multiplicity distributions
From the identified protons and anti-protons, event-by-event net-proton distribution are constructed. Figure 14 Fig. 15 is found to be smaller than that observed for 0-5% centrality.

Collision energy dependence of net-proton cumulants
Cumulants up to the 4 th order of the event-by-event net-proton multiplicity distributions for most central (0-5%) Au+Au collisions as a function of collision energy ( √ s N N ) are shown in Fig. 16. The data  points in blue are the measurements corrected for detector efficiency and acceptance effects while those in red, represent the efficiency uncorrected results. All the presented results are corrected for volume fluctuation effects by employing the CBWC procedure. The efficiency corrected cumulant C 1 and C 3 shows monotonic decrease as a function of collision energy, whereas a non-monotonic collision energy dependence is observed for efficiency corrected C 2 and C 4 . The trend of collision energy dependence shown by the efficiency corrected cumulants is already captured by the efficiency uncorrected cumulants shown in the open red marker. Upto third order of cumulant (C n , n ≤ 3), the efficiency uncorrected cumulants remain consistently lower than efficiency corrected results at all energies. No such hierarchy is observed for C 4 . The efficiency correction increases the magnitude of statistical uncertainties on the cumulants. Figure 17 shows the collision energy dependence of efficiency corrected and uncorrected net-proton cumulants for the peripheral 70-80% collisions. The efficiency corrected cumulants up to the fourth order is larger than their corresponding efficiency uncorrected counterparts across all the energies. The collision energy dependence trend of cumulants for peripheral 70-80% collisions is weaker as compared to the 0-5% central collisions. As seen from the peripheral measurements, the magnitude and collision energy dependence trend of odd order cumulants (C 1 and C 3 ) are quite similar. The same observation can also be made for C 2 and C 4 . This is a characteristic of a Skellam distribution where all the odd order cumulants are the same and are simply the difference of the mean of constituent Poisson distribution, and all the even order cumulants are also the same and given by the sum of the means.

Collision energy dependence of net-proton cumulant ratio
Recent measurements on net-proton cumulant ratios up to fourth -order have shown an intriguing collision energy dependence trend for most central 0-5% Au+Au collisions. The collision energy dependence of higher-order cumulant ratio C 4 /C 2 , will be discussed in this subsection in particular. As discussed  previously, to eliminate the trivial system volume dependence to allow a direct comparison to theoretical calculations, ratios of cumulants are constructed. Figure 18, shows the collision energy dependence of cumulant ratio C 4 /C 2 (equivalently moment product κσ 2 ) of net-proton distributions in central 0-5% and peripheral 70-80% Au+Au collisions at RHIC. Various variants of HRG and UrQMD model calculations are also shown as baselines to the measurements. The STAR net-proton κσ 2 measurements for the most central 0-5% collisions show a non-monotonic collision energy dependence, which is neither reproduced by the transport model UrQMD nor the various variants of HRG model calculations presented. The ideal-HRG (HRG GCE in the figure), remains close to the Skellam baseline value at unity. The HRG calculations in the grand canonical framework taking only repulsive interactions into account are given by the red dotted line. The calculations remain below unity, and the deviation from unity increases with decreasing collision energy. With both repulsive and attractive interactions taken in HRG calculations, the κσ 2 shows a larger deviation from unity towards lower values as seen from the black dotted line. The UrQMD calculation for the most central 0-5% collisions suggests a monotonic decrease of κσ 2 with decreasing collision energy. This suppression of cumulant ratio C 4 /C 2 in the lower collision energies could be attributed to the effect of baryon number conservation and is also reproduced by HRG canonical ensemble calculations [83]. Towards the higher collision energies, measurements show agreement with lattice-QCD calculations in the energy range √ s N N = 39 -200 GeV [65] subjected to caveats discussed in sub-section 3.1. Interestingly, FRG calculation [61] shows a decreasing trend from higher to lower energy up to √ s N N = 10 GeV, below which it shows an increasing trend, which the authors [61] attribute to the sharpening of crossover. κσ 2 from peripheral collisions, on the other hand, shows a weak collision energy dependence and always remains smaller than the Skellam baseline at unity. The observation of non-monotonic collision energy dependence of κσ 2 for most central 0-5% collisions is qualitatively consistent with expectations from the linear sigma model calculations, which  Fig. 19 as a function of collision energy. Various baselines, starting from the Skellam baseline at unity to UrQMD and HRG model calculations, are considered. Also shown is the case where the measurements at peripheral 70-80% collisions serves as a baseline to measurements at central collisions. Irrespective of which baseline one chooses, the most central 0-5% measurements show both positive and negative deviations from the baseline across collision energies at a level of 2-3σ. This suggests the robustness of the non-monotonic collision energy dependence with respect to various baseline calculations. This observation is further strengthened by performing a χ 2 test to demonstrate the deviation of the experimental data to various model calculations in the collision energy range √ s N N = 7.7 -27 GeV [31,32]. The difference between experimental data and model expectations is quantified in terms of χ 2 and then converted to right-tailed p-value. The right-tailed p-value is obtained as: p-value = 1 − cdf χ 2 ,d (χ 2 c ), where cdf χ 2 ,d denotes the cumulative distribution function of the χ 2 distribution with d degrees of freedom and χ 2 c is the calculated χ 2 value. From the table of the probability distribution of χ 2 [119], right-tailed p-value can be obtained for a given number of degrees of freedom. A criterion on the right-tailed p-value: p < 0.05 is commonly used to ascertain a significant deviation between the data and model expectations. The p-value from the χ 2 test between data and various variants of HRG and UrQMD model calculations are tabulated in table 3. The small p-values (p < 0.05) obtained from the test demonstrate that the data significantly differ from model calculations.
The observed non-monotonicity has been quantified by using a statistical procedure. The most central 0-5% measurements on net-proton κσ 2 as a function of collision energy, in the collision energy range √ s N N = 7.7 -62.4 GeV, is best fitted by a polynomial of degree four with a χ 2 per degree   GeV is obtained to be 3.1σ based on the fourth order polynomial fitting procedure [31]. If the baseline is changed to UrQMD and HRG canonical ensemble calculation, the fourth order polynomial fit to data results in χ 2 per degree of freedom ∼1.0. The sign change of derivative for the UrQMD and HRG CE baselines are shown in the bottom panels, and using the statistical procedure as used in Ref. [31], the significance of non-monotonic dependence of measurements with respect to both these baselines are found to be 3.3σ.
Measuring the proton and anti-proton cumulants could reveal what drives the observed non-monotonic  The measurements are shown for 0-5% and 0-80% collision centrality class. The bars and brackets on the data points represent the statistical uncertainties and systematic uncertainties, respectively. UrQMD and JAM model expectation [118] for 0-5% central Au+Au collisions are given by the green and brown bands, respectively. Various variants of HRG model calculations are also presented. Lattice-QCD and FRG calculations are shown as red and blue bands, respectively [61,65].
collision energy dependence of net-proton cumulants. The cumulant ratios of proton and anti-proton distributions subtracted from unity as a function of collision energy for most central (0-5%) and peripheral (70-80%) Au+Au collisions at RHIC are given in Fig. 21. The second and third-order cumulant ratios minus unity remain negative across all collision energies for both protons and anti-protons. For protons, the deviation of the measurements from Poisson baseline at zero increases with decreasing collision energies, whereas for anti-protons, the deviations from Poisson baseline is maximum at √ s N N = 200 GeV. This is understood to be coming from two-particle correlation functions of protons and anti-protons [32]. The qualitative trend of both ratios is captured by the UrQMD model for protons and anti-protons. C 4 /C 1 − 1 of proton distributions shows a non-monotonic collision energy dependence similar to that observed in net-proton cumulant ratio C 4 /C 2 , whereas C 4 /C 1 − 1 of anti-proton distributions show weak dependence on collision energy. The UrQMD calculations show a monotonic collision energy dependence of C 4 /C 1 − 1 for protons, whereas for anti-protons, it is mostly consistent with Skellam baseline at zero. It was found in ref [32], that the non-monotonic collision energy dependence of proton C 4 /C 1 − 1 is driven by the presence of four-particle correlations for protons, whereas for anti-protons, only the two-particle correlations drive the collision energy dependence trend.

3.4.4
Collision energy dependence of net-charge and net-kaon cumulant ratio C 4 /C 2 The ratio of fourth to second-order cumulant of net-charge and net-kaon distributions in most central (0-5%) and peripheral (70-80%) Au+Au collisions at STAR are shown in Fig. 22. For net-charge cumulant  measurements, all charged particles within p T range 0.2 -2.0 GeV/c and |η| < 0.5 are selected and the centrality is determined uisng charged particle multiplicity in the region 0.5 < |η| < 1.0. Net-kaon cumulant measurements are done selecting charged kaons (K + and K − ) within 0.2 < p T < 1.6 GeV/c at mid-rapidity (|y| < 0.5). All charged particles in the |η| < 1.0 coverage excluding kaons (K + and K − ) are used in the centrality definition. The most central (0-5%) net-charge and net-kaon cumulant ratio C 4 /C 2 have flat collision energy dependence. The measurements, inlcuding those from the PHENIX experiment, are mostly consistent with the Skellam baseline at unity. The statistical uncertainties on the measurements are large, especially at the lower collision energies except for the measurement from the PHENIX experiment, where they are small due to smaller detector acceptance [120]. For the net-charge fluctuation measurements, the large statistical uncertainties are attributed to the larger width of net-charge distributions as compared to the other two net-particle measurements, whereas smaller detection efficiency of charged kaons contributes to large statistical uncertainties on net-kaon cumulant ratio C 4 /C 2 . The larger width of the net-charge distributions can be attributed to the decay of resonances. The net-charge C 4 /C 2 in peripheral (70-80%) collisions remains larger than the Skellam baseline across all the energies, whereas the net-kaon C 4 /C 2 for the same collision centrality are either consistent or higher than Skellam baseline at unity, at all the energies. One should also keep in mind that in addition to charged kaons, neutral kaons and multi-strange baryons like the Λ are also produced in appreciable numbers in heavy-ion collisions and carry strangeness. Thus, more strange particles should be added to net-kaon fluctuations to serve as a more effective proxy of net-strangeness fluctuations. The lack of statistical precision in these measurements from BES-I forbids us to see any possible critical signal at present. The second phase of the BES program, BES-II, with large event statistics, will shed more light on this.

Acceptance dependence of net-proton cumulant ratios
The acceptance dependence of number fluctuations are realized by changing the pseudo-rapidity (η) or rapidity (y) coverage considered for fluctuation measurements. Such dependence could reveal important information on the nature of the fluctuations and their origin. The fluctuations could be broadly divided into two categories depending on their source of origin: (a) Thermal fluctuations and (b) Initial state fluctuation. While the thermal fluctuations increase with acceptance and finally saturate after the typical equilibrium correlation length (ξ ∼ 0.5f m) is reached, the initial state fluctuations are longrange in nature and grow with acceptance till the effect of conservation laws become significant. In a very small acceptance window, the fluctuations approach the Poisson limit while at broader rapidity coverage, they are driven by conservation laws. Also, it has been suggested in Ref. [121], rapidity dependence of fluctuations could also help understand the diffusion property of the system formed in the heavy-ion collisions. If the critical point lies between two energies in the beam energy scan, then studying the rapidity dependence of cumulants can help pinpoint its location [122]. Rapidity dependence of cumulants measured by the STAR and ALICE experiments with the limitations of finite detector acceptance are discussed in this section. The rapidity dependence of net-proton C 4 /C 2 in the collision  Fig. 23 [31]. The ratio remains close to the Poissonian limit at unity for the smallest rapidity window accessed at all energies. While decreasing trend of C 4 /C 2 with increasing rapidity acceptance is observed for 19.6 and 27 GeV, the ratio at 7.7 GeV first shows a dip followed by an increasing trend, albeit with large uncertainties when acceptance is enlarged. For the higher collision energies √ s N N ≥ 39 -200 GeV, weak dependence is observed on acceptance. ALICE also reported its first results on net-proton cumulant measurements up to the second order. The ratio of second order net-proton cumulant to that of mean value of protons+anti-protons are shown as a function of η coverage in Pb+Pb collisions at √ s N N = 2.76 TeV [123]. Weak acceptance dependence of measurements is observed for second-order ratio, with central values of the measurements differing at a level of 4% from lowest to largest acceptance. The measurements approach the Poisson limit when acceptance is reduced. The HIJING model calculations fails to describe the measurements. The data seem to be consistent with the global baryon number conservation scenario rather than local baryon number conservation. Increasing the p T range of the charged particle selection in the measurements could also provide interesting results as more charged particles now constitute the system. The recent result on non-  monotonic energy dependence of net-proton C 4 /C 2 [31] was realised after the p T range for protons was enlarged to 0.4 < p T < 2.0 GeV/c as compared to the prior measurements with 0.4 < p T < 0.8  GeV/c [30]. Extending the p T range was possible due to the addition of the TOF detector. Figure 25 shows the p T dependence of net-proton C 4 /C 2 in most central 0-5% Au+Au collisions at √ s N N = 27 GeV. The lower limit of p T is fixed at 0.4 GeV/c so as to avoid contributions from background protons from beam pipe interactions, and upper limit is varied. The data points show a larger deviation from the Skellam baseline at unity with increasing p T coverage, and the statistical uncertainties on the data points increase as more particles are selected with larger p T acceptance. The critical phenomenon are of long-wavelength character and thus should be predominant in the lower momenta. Reducing the current lower limit of p T = 0.4 GeV/c is a highly desired and challenging task in the search of the CP. Increasing the phase space coverage by reducing the lower limit of p T would also increase the possibility of ∆ and weak decays of strange baryon induced self correlations, as discussed in subsection 3.3.6. By performing a similar study as shown in Fig. 13, but systematically reducing the lower limit of p T down to 0.05 GeV/c, a negligible self-correlation was found due to such decays.

A Lattice-QCD inspired fit to experimental measurements
In the previous subsection, the cumulant ratios measured at STAR were compared with various model calculation. This section discusses the comparison of net-proton cumulant ratios C 3 /C 2 to a fit inspired by lattice-QCD calculations for net-baryon number fluctuations. The net-proton cumulant ratios Sσ = C 3 /C 2 from nine collision energies as a function of M/σ 2 = C 1 /C 2 is presented in Fig. 26. As discussed in Ref. [19], in lattice-QCD calculations, the baryon number susceptibility ratios χ B 3 /χ B 1 and χ B 3 /χ B 2 can be expressed in terms of χ B 1 /χ B 2 as follows.  The quantities p a1 and p b1 could be fixed by comparing the experimental data on net-proton C 3 /C 1 as a function of C 1 /C 2 to eqn. 71 for the collision energy range √ s N N = 7.7 -200 GeV [31] assuming net-proton to be proxy for net-baryon. The p a1 and p b1 are then used to calculate C 3 /C 2 as a function of C 2 /C 1 using the relation mentioned in eqn. 72. The predicted dependence is shown as green band in Fig. 26. The dashed line at y=x refers to the case where the ratio Sσ = M/σ 2 . The HRG model calculations is found to be along the dashed line whereas the UrQMD expectations show small deviation at large M/σ 2 corresponding to low collision energies. The STAR measurements lie on the right side of the dashed line at y=x, indicating that M/σ 2 > Sσ at all energies. This hierarchy is consistent with LQCD calculations reported in Ref. [65]. The agreement of the experimental measurements to the LQCD inspired fit within uncertainties suggests the production of strongly interacting and thermalised QCD matter in heavy-ion collisions.

3.4.7
Other new measurements   [124,65]. The UrQMD calculations are either positive or consistent with zero. Collision energy dependence of fifth and sixth-order fluctuations has also been reported by STAR [125]. Collision energy dependence of netproton cumulant ratio C 5 /C 1 and C 6 /C 2 is shown in Fig. 28 and 29, respectively [125]. LQCD, FRG calculations along with expectations from the UrQMD model and HRG CE model are also presented. To improve the statistical precision of the measurements, calculations are done over 0-40% centrality. The results for most peripheral 70-80% centrality are also presented. The C 5 /C 1 shows weak collision energy dependence. Deviations of the measurements from zero at a level of ≤ 2σ significance are observed. The LQCD and the FRG calculations show negative values of the ratio of fifth-to-first order baryon number susceptibilities. The peripheral 70-80% measurements, on the other hand, remains consistently positive across all collision energies and close to the Skellam baseline at unity. The net-proton C 6 /C 2 in central 0-40% collision shows increasing negative values with decreasing collision energy. This is qualitatively consistent with the sign and trend predicted by LQCD and FRG calculations for the ratio of sixth-tosecond order baryon number susceptibilities. The magnitude of deviation of the measurements from zero is of similar significance, as for C 5 /C 1 . In contrast, the peripheral 70-80% centrality measurements are positive at all the energies (shown in the inset of Fig. 29). The HRG CE expectation for both the net-proton cumulant ratios are positive in the higher collision energies and show suppression towards negative values in the lower collision energies due to the effect of baryon number conservation. Furthermore, a particular ordering of net-proton cumulant ratios: C 3 /C 1 > C 4 /C 2 > C 5 /C 1 > C 6 /C 2 as predicted by lattice-QCD was also observed in data for Au+Au collisions at √ s N N = 7.7 GeV and 200 GeV. The ordering of the four ratios C 3 /C 1 , C 4 /C 2 , C 5 /C 1 and C 6 /C 2 for net-proton distributions in 0-40% collision centrality are shown in the Fig. 30. While the FRG calculations also follow the predicted hierarchy, the UrQMD does not show such ordering of ratios and gives similar positive values for the four ratios. Here, one should keep in mind the caveats when comparing experimental data to LQCD and FRG model calculations. As opposed to the LQCD and FRG, which calculates net-baryon fluctuations, the experimental measurements are made for net-proton fluctuations within a specific kinematic phase space allowed by the experiment. Also, the former does not take into account the dynamics associated with nuclear collisions, for example, the collective expansion. These limitations should be accounted for when making a quantitative comparison between data and theory.   [65], FRG [61] and UrQMD calculations (0-40% centrality) are shown as red, grey and brown bands, respectively.
The fifth and sixth order factorial cumulants (κ 5 and κ 6 , respectively) of proton multiplicity distributions as a function of collision centrality in Au+Au collisions at √ s N N = 7.7 GeV are presented in Fig. 31 [125]. The κ 5 shows a decreasing trend with increasing collision centrality, whereas a weak collision-centrality dependence is found for κ 6 . The UrQMD model calculations for both the measurements are close to the Poisson baseline at zero. The κ 5 measurement for 0-5% centrality is consistent with the expectation from a two-component proton multiplicity distribution model, which takes proton cumulants up to fourth-order as input in its construction. The most central (0-5%) κ 6 , however, is 1.8σ away from expectation from such a model and consistent with the Poisson baseline at zero, albeit with large uncertainties. These higher-order cumulant measurements are expected to be measured with higher precision in the second phase of the beam energy scan program. (b) Multiplicity dependence of net-proton cumulants in pp collisions : Cumulant measurements have also been recently explored in small systems by STAR experiment [126]. Multiplicity dependence of net-proton cumulant ratios C 4 /C 2 , C 5 /C 1 and C 6 /C 2 in p+p collisions at √ s = 200 GeV are shown in Fig. 32. The three ratios approach to their Skellam and HRG GCE baseline value of unity with increasing charged particle multiplicity while the centrality dependence measurements of these ratios in Au+Au collisions at √ s N N = 200 GeV show opposite trend. CBWC average of p+p collision results are consistent with the multiplicity dependence trend. The C 5 /C 1 and C 6 /C 2 in p+p collisions are positive across the entire charged particle multiplicity range which is in contrast to the negative sign predicted by LQCD calculations for net-baryon χ B 5 /χ B 1 , χ B 6 /χ B 2 for QCD matter. On the other hand, the C 5 /C 1 and C 6 /C 2 for Au+Au collisions at √ s N N = 200 GeV show increasingly negative value towards central collisions. The Pythia calculations are found to overestimate the values of the ratios.
The higher-order cumulant measurements in a small system provide a baseline for comparison with results from heavy-ion collisions. They also provide important data to test thermalization and freeze-out conditions.
(c) Net-Λ fluctuations: Net-Λ fluctuations have been studied at RHIC in the energy region √ s N N = 19.6 -200 GeV [127]. The collision energy dependence of net-Λ fluctuations is shown in Fig. 33. Monotonic collision energy dependence is observed for the cumulant ratios C 2 /C 1 and C 3 /C 2 in both central (0-5%) and peripheral (50-60%) collisions. NBD baseline seems to show better agreement with  Comparison of net-proton measurements with those of net-Λ suggests that the net-Λ results follow the qualitative energy dependence as shown by net-proton. Experimentally Λ particles are detected via invariant mass reconstruction of their decay daughters, protons, and pions. To treat the net-Λ and net-proton fluctuations as independent, experimental measurements should ensure the decay contributions from Λ do not contaminate the net-proton sample. If this can be achieved either by applying kinematic cuts for particle selection or feed-down corrections, then these two net-particle fluctuations could be added to construct a more reliable proxy of net-baryon fluctuations. Since, Λ has both baryon and strangeness quantum numbers, their inclusion is essential in the baryon-strangeness correlation studies. Studies on the non-diagonal baryon-strangeness correlator by authors of Ref. [128], indicate that considering only net-protons and net-kaons in such correlation studies are not sufficient. In fact, proton-kaon correlation has negligible contribution in the baryon-strangeness correlations, and a significant contribution to the latter comes from the variance of net-lambda distribution. In future, it will be very interesting to see if the C 4 /C 2 for net-Λ as a function of collision energy shows a similar non-monotonic variation as observed for net-protons.
(d) Deuteron and nuclei fluctuations : Since deuteron is the second lightest stable baryon detected in the heavy-ion collision experiments, measurement of their fluctuations might have interesting prospects in the search for the QCD critical point. Especially, in the high baryonic density region (low collision energies), an abundance of the light nuclei is enhanced, and their contribution to the dynamics of baryonic charge fluctuations in the system becomes important. Higher-order moments of event-by-event number distribution of deuterons are recently suggested also as a probe for their production mechanism [129]. Synthesis mechanism of light nuclei detected in high-energy nuclear collisions is  primarily discussed in two scenarios: (a) statistical thermal model and (b) coalescence model. Though both scenarios, to a good extent are able to explain the experimental data on yields and their ratios, a complete picture of production mechanism of light nuclei is still missing. In this context, measurement of light nuclei fluctuations is worthwhile to probe their production mechanism and explore signals for CP [130].
(e) Probing magnetic field via fluctuations : Presence of a strong magnetic field in heavy-ion collisions can bring interesting effects on QCD phase structure. The chiral-crossover transition temperature T pc is expected to be reduced due to the phenomenon of inverse magnetic catalysis [131,132,133,134]. Lattice-QCD calculation studies done by authors of Ref. [135] suggest using fluctuations and correlation of conserved charge quantities to probe the existence of magnetic field in the heavy-ion collisions. It has been known from previous studies, in the absence of a magnetic field, the second-order fluctuations of B, Q, and S increase with increasing temperature [136,65]. Recent studies indicate that they seem to develop a peak structure due to the presence of a strong magnetic field, with the location of the peak decreasing towards lower temperatures with increasing magnetic field. In particular, observables constructed from baryon-strangeness, charge-strangeness and baryon-charge second order correlations and fluctuations: and −3χ BS 11 /χ S 2 have been proposed to detect presence of magnetic field in early stages of heavy-ion collision and breaking of isospin symmetry in non-zero magnetic field [135]. The proposed magnetic-field-dependence of the observables could be realized by performing measurements on the centrality dependence of these quantities in experiments as the magnetic field is expected to vary as a function of collision-centrality.

Future directions
Discovering the QCD critical point will be a landmark in the QCD phase diagram. Progress in both theory and experimental measurements is needed to ascertain the existence and location of the QCD critical point. We discuss in this section the outlook from the theoretical side on the study of critical fluctuations along with the future opportunities with upcoming new experimental facilities in the hunt for the QCD critical point.

Theory
A concrete answer to the existence and location of the QCD critical point will require a complete understanding of the dynamics of the system created in the heavy-ion collisions. A quantitative theoretical framework that concerns the dynamical evolution of the system formed in heavy-ion collisions till it reaches its final state will be crucial to the search of QCD critical point. Studies of fluctuations and their coupling with critical modes in stochastic-hydrodynamics as well as hydro-kinetics have been ongoing. Stochastic fluid dynamics generally involves event-by-event simulation of the viscous relativistic fluid dynamical evolution with stochastic conservation law, where energy-momentum tensor and conserved current 4-vector include Ideal + Viscous + Noise parts. On the other hand, the hydro-kinetics method is used to calculate the time dependence of the two-point correlation functions of the fluid dynamical fields using a set of linearized equations of motion. Hydrodynamical fluctuations and their time evolution considering the effect of critical slowing down are being studied [138,139,140,141,142,143,144,145]. Another important aspect of stochastic fluid dynamics and hydro-kinetics study of fluctuations is the particlization. Experimental observables measure correlations of produced particles. The conversion from fluid degrees of freedom to the particle degrees of freedom includes the effects like finite particle statistics, limited acceptance and conservation of energy and global charges at the freeze-out hypersurface.
Further, the thermodynamics of the system formed in these collisions can be studied via comparison of the fluctuations measured in the experiments to the various thermal model and lattice-QCD calculations [26,146,147]. Since the net-particle distribution could be fully characterized by its all-order moments, consistent thermal description of mean yields as well higher moments becomes important to test the thermal nature of the whole system [148]. Deviation from thermal equilibrium might indicate interesting physics, such as the presence of a QCD critical point.

Experiment
The observed non-monotonic collision energy dependence of net-proton C 4 /C 2 has been the most interesting result so far in the experimental search of the QCD critical point. Recently STAR collaboration reported the proton C 4 /C 2 in Au+Au collisions at √ s N N = 3 GeV from STAR's FXT program [149]. At such low energies, the number of event-by-event anti-protons are negligibly. As seen from the Fig. 34, a suppression in C 4 /C 2 was observed compared to higher collision energies. This is consistent with expectations from baryon number conservation and suggests that matter is dominantly hadronic at such low collision energies. The non-monotonic collision energy dependence of the C 4 /C 2 measurements from the BES-I program together with the observed suppression at √ s N N = 3 GeV (FXT) indicates that the possible critical point could only exist at collision energies higher than 3 GeV. Nonetheless, the measurements at collision energies √ s N N = 7.7 -27 GeV suffer from large statistical uncertainties which are the key energies that drive the non-monotonic collision energy dependence trend. The measurement on proton C 4 /C 2 by HADES [105], for 0-10% central Au+Au collisions at √ s N N = 2.4 also has large uncertainties and require more event statistics to serve as effective experimental baseline in the large µ B region. Here, in this section, we discuss future prospects of fluctuation measurements in experiments that are currently active and upcoming in the near future.

STAR BES-II and FXT
In order to make precise measurements in the collision energy range 7.7 GeV ≤ √ s N N ≤ 19.6 GeV, phase II of the beam energy scan (BES-II) has been ongoing using the STAR detector at RHIC [150]. As compared to the presented results from the BES-I program, additional collision energy at √ s N N = 9.2 GeV from the BES-II program will be available for fluctuation measurements with large statistics. Table 4 contains the event statistics for each collision energy being collected in the ongoing BES-II program and the corresponding µ B at chemical freeze-out. The statistical uncertainties of the measurements in BES-II will be significantly reduced as compared to current measurements. The magnitude of projected statistical uncertainties from BES-II measurements are shown in Fig. 34. In addition to precision measurements of cumulants, the goal is to also improve the quality of the potential critical signal. Several detector upgrades have been made in the STAR detector at RHIC. Upgrades such as that of the inner chambers of the TPC, the event plane detector, the end-cap time of flight detector, bring new scopes to the fluctuation measurements, like enlarged rapidity coverage (from |η| < 1 to |η| < 1.5), improved particle identification and centrality definition at forward rapidity. The QCD critical point induced power law behaviour of the fluctuation measurements can be tested using the extended rapidity coverage [151]. STAR fixed target (FXT) program will carry Au+Au collisions from √ s N N = 3.0 -7.7 GeV (shown in the Fig. 34 [154], and (d) Heavy-Ion program at Japan Proton Accelerator Research Complex (JPARC-HI), Tokai, Japan [155].

NICA
Two dedicated experiments are designed at NICA for the study of QCD phase diagram: a fixed target experiment, Baryonic Matter at Nuclotron (BM@N) [156], which will cover the lower energy region, and the Multi-Purpose Detector (MPD) [157] which will run in collider mode. The MPD experiment includes critical point search using fluctuation measurements among its physics priorities [158]. The MPD is a 4π spectrometer which can detect charged hadrons, electrons and photons in heavy-ion collisions in the energy range of the NICA collider. The detector setup includes the Central Detector (CD) covering ±2 units in pseudo-rapidity (η), Time-Projection Chamber (TPC), barrel Time-Of-Flight system (TOF), Zero-Degree Calorimeter (ZDC) and Fast Forward Detector (FFD) [159]. The TPC is the main tracking detector of the MPD and will provide tracking and PID with resolution better than 8% in the range |η| < 1.2, and momentum resolution for charge particles better than 3% in 0.1 < p T < 1 GeV/c. The cylindrical part of TOF also has full azimuthal and |η| < 1.2 coverage and will allow PID in the p T range 0.1 -2 GeV/c for charged hadrons. Combining the measurements from TPC and TOF will provide an efficient π/K PID separation up to 1.5 GeV/c and π/p separation up to 3 GeV/c [160]. Initial run of NICA will be done with Bi+Bi collisions at √ s N N = 9.2 AGeV [161]. Au+Au collisions are planned with beam energies spanning from √ s N N = 12 -27 GeV with luminosity L ≥ 10 30 cm −2 s −1 for minimum bias is planned for the higher collision energies but decreases to about 10 Hz because of low luminosity at √ s N N = 4 GeV. A technique to sustain the beam luminosity in low energies for prolonged period of time, called the electron cooling, which is also being used in BES-II at STAR, will be employed.

Compressed Baryonic Matter (CBM) experiment
The CBM experiment is a fixed-target detector capable of identifying hadrons, leptons (e, µ) and photons. The interaction rate will be higher for the CBM experiment at FAIR. The experiment will operate over the energy range, √ s N N = 2.7 -4.9 GeV. The goal of the CBM experiment at SIS100 is to discover fundamental properties of QCD matter at high baryon density using rare probes. The experiment plans to perform high-precision study of higher-order fluctuations at various beam energies in baryon density region of 500 MeV < µ B < 800 MeV. In order to achieve high precision, the measurements will be performed at reaction rates up to 10 MHz. The CBM detector acceptance of polar emission angles between 2.5 and 25 degrees allows it to cover mid-rapidity and the forward rapidity for symmetric collision systems over the FAIR energy range. The lab pseudo-rapidity range coverage of the CBM detector system is 1.5 < η < 3.8. The detector setup includes Micro Vertex Detector (MVD), Silicon Tracking System (STS), Ring Imaging Cherenkov Detector (RICH), Transition Radiation Detector (TRD), and Time-of-Flight Detector (TOF) [162,163]. The collision vertex of the event will be provided by the MVD. The STS inside the dipole field records the momentum information of the charged particle tracks. Particle identification will be done by the TOF. Measurement of centrality of collision will be possible with Project Spectator Detector (PSD). Feasibility study of performing fluctuation measurements at CBM are ongoing and show promising prospects [164].
The new experiments will collect enormous event statistics due to high interaction rates. Large event statistics will be beneficial for higher-order fluctuation measurements. Handling such rates would demand fast and state-of-the-art detectors, which are currently being envisaged. The interaction rates from various upcoming experiments are summarised in the Fig. 35 [165,55]. In addition, understanding the systematics of fluctuation measurements is crucial. For example, determining the centrality from different detectors (PSD at CBM and EPD at STAR) than those used for PID, enlarging the kinematic acceptance for PID (iTPC at STAR) will help us to study the measurements better. Due to low multiplicity of produced particles in low energy collisions, centrality definition using charge particles or energy deposition measurements suffer from poor resolution. This could result in a artificial correlation in fluctuation measurements as discussed in subsection 3.3.5. So, proper investigation of this effect should be studied in the fluctuation measurements. As nearly half of the net-baryon is carried out by (anti-)neutrons, fluctuations of neutron number, in principle, should also reflect criticality. Adding neutrons to proton fluctuations would serve as a better proxy of net-baryon number. Other important factors will be the detector efficiencies and acceptance for various hadrons at low energies, especially in the fixed target configuration. Uniform acceptance of protons at mid-rapidity coverage would be crucial for studies on the critical point via proton number fluctuations. In nuclear collision experiments, for a given detector setup, the mode of collisions determines the phase space acceptance for the charged particles. For example, as shown in Fig. 36, the acceptance of protons is asymmetric in fixed target mode (upper panel) as compared to the collider mode (lower panel), where it is symmetric.

Acceptance: collider vs fixed target
In the collider mode, the mid-rapidity is always at zero, and this allows STAR TPC to perform mid-rapidity measurements at all collision energies. However, in the fixed target mode, the location of mid-rapidity depends on collision energy. The acceptance of the STAR detector allows coverage of mid-rapidity distributions of protons in fixed target mode at √ s N N = 3 GeV, while at higher fixed target energies, additional detectors like iTPC (inner -TPC) and eTOF (extended -TOF) are needed for the phase-space coverage around the mid-rapidity.

Summary
In the absence of any first-principle theoretical calculations on the location of the QCD critical point, one has to rely on predictions from QCD-based theories to search for critical signals in heavy-ion collisions. One such calculation is done with the linear sigma model, which predicts a non-monotonic dependence of kurtosis of number fluctuations when approaching the critical point. Experimental measurement on the net-proton cumulant ratio C 4 /C 2 in 0-5% Au+Au collisions, shows an evidence of a non-monotonic collision energy dependence of the measurements in the range √ s N N = 7.7 -200 GeV.
The normalised fourth-order cumulant ratio of proton shows collision energy dependence similar to net-proton measurements, which is attributed to the presence of four-proton correlations, the antiproton, on the other hand, have a weak dependence on collision energy. None of the model calculations are able to reproduce the observed energy dependence of the net-proton cumulant C 4 /C 2 . Comparison with a wide pool of baseline calculations suggests that the 0-5% net-proton C 4 /C 2 measurements deviate below and above the baseline calculations at a level of 2-3 σ at various collision energies. The exhibited non-monotonicity of net-proton cumulants with respect to Skellam baseline, UrQMD baseline, and HRG canonical ensemble calculations are found to be of 3.1-3.3 σ significance. The recent measurement of proton C 4 /C 2 in Au+Au collisions at √ s N N = 3 GeV show suppression for 0-5% centrality and is consistent with the expectation from the baryon number conservation effect suggesting the dominance of hadronic matter in such low energies. The suppression of proton cumulant ratio C 4 /C 2 at 3 GeV, together with the reported non-monotonic collision energy dependence of net-proton C 4 /C 2 at higher energies, indicates that the possible critical point could exist at collision energies higher than 3 GeV. The experimental results on net-charge and net-kaon cumulant ratio C 4 /C 2 in most central 0-5% collisions have large uncertainties associated with the measurements and show a flat collision energy dependence with most of the data points showing agreement with the Skellam baseline within uncertainties. The current net-Λ fluctuation results are only up to the 3 rd order and do not shown any signs of critical fluctuations. Addition of net-Λ's to net-kaon fluctuation could serve as a better proxy for net-strangeness fluctuations as most of the strangeness is carried by kaons and Λ's in the system.
Negative values were observed from measurements of the fifth and sixth-order net-proton cumulant ratios, C 5 /C 1 and C 6 /C 2 for 0-40% Au+Au collisions at most of the collision energies, although with large uncertainties. C 6 /C 2 becomes increasingly negative with decreasing collision energy. At √ s N N = 200 GeV, the C 6 /C 2 measurements progressively become negative from peripheral to central collisions. These observations, albeit with large uncertainties, are qualitatively consistent with the trend and sign predicted by LQCD calculations for the ratios of fifth-to-first and sixth-to-second order baryon number susceptibilities for QCD matter in the range µ B < 110 MeV. In contrast, positive C 5 /C 1 and C 6 /C 2 were reported in peripheral (70-80%) collisions at all energies. Multiplicity dependence studies of the two ratios in p+p collisions at √ s N N = 200 GeV also yeild positive sign. Current measurements on higher order factorial cumulants at √ s N N = 7.7 GeV indicate that more event statistics is needed for the search of first-order phase transition in the high µ B (∼ 400 MeV) region. The system formed in the heavy-ion collision is of femto-scale level and evolves dynamically with time. Hence, it is quite challenging both theoretically and experimentally to probe into the physics of QCD critical point. Theoretical calculations considering the effects of non-equilibrium and critical slowing down near the CP will guide the experimental search for the CP. As a first step, one tries to establish the possibility of the existence of the QCD critical region, which could span about 100 MeV range in baryon chemical potential. The current result on net-proton kurtosis is an important step in that direction as these are the first fluctuation measurements exhibiting non-monotonic variation with collision energy. To reduce the statistical uncertainties on the current net-particle cumulant ratios measurements at lower collision energies and improve the quality of the observed signal, the BES-II program has been carried out at RHIC. The STAR FXT program, along with several future experiments like, NICA and CBM, will collect large event statistics and explore the baryon-rich energy region. The overlapping µ B region covered by these experiments will allow for independent confirmation of possible critical signals. With the state-of-the-art detector technology at the helm, the upcoming experiments carry huge potentials to discover the elusive critical point and map the QCD phase diagram with utmost precision.