# Correlated fluctuations near the QCD critical point

## Abstract

In this paper, we introduce a freeze-out scheme for the dynamical models near the QCD critical point through coupling the decoupled classical particles with the order parameter field. With a modified distribution function that satisfies specific static fluctuations, we calculate the correlated fluctuations of net protons on the hydrodynamic freeze-out surface. A comparison with recent STAR data shows that our model calculations could roughly reproduce energy dependent cumulant and of net protons through tuning the related parameters. However, the calculated and with both Poisson and Binomial baselines are always above the experimental data due to the positive contributions from the static critical fluctuations. In order to qualitatively and quantitatively describe all the related experimental data, the dynamical critical fluctuations and more realistic non-critical fluctuation baselines should be investigated in the near future.

###### pacs:

05.70.Jk, 25.75.Gz, 25.75.-q, 25.75.Nq## I Introduction

.

The Beam Energy Scan (BES) program at the Relativistic Heavy Ion Collisions (RHIC) aims to explore the QCD phase diagram and study the QCD phase transition Aggarwal:2010cw (). At small chemical potential region, lattice simulations demonstrate that the QCD phase transition is a rapid but smooth crossover Karsch:2003jg (); Aoki:2006we (); Aoki:2009sc (); Bazavov:2014pvz (). At large chemical potential and medium temperature region, different effective models show that the QCD phase transition is the first order with a clear phase boundary Klevansky:1992qe (); Fukushima:2003fw (); Roberts:1994dr (); Berges:2000ew (). Correspondingly, a critical point is located at the end of the first order phase boundary, leading to unique critical phenomena of the QCD matter Stephanov:2004wx ().

In the vicinity of the critical point, the correlation length and the magnitude of fluctuations become divergent in a static and infinite medium. Detailed calculations have shown that the correlated fluctuations between two particles are proportional to the square of the correlation length: Stephanov:1998dy (); Stephanov:1999zu (), higher cumulants of the correlated fluctuations are more sensitive to the correlation length with and Stephanov:2008qz (); Athanasiou:2010kw (). It was also found that the fluctuations of net protons are related to the baryon susceptibilities Hatta:2003wn (); Kitazawa:2012at () - the quantities that are generally used to evaluate the critical fluctuations in Lattice QCD and effective field theories. Therefore, the higher cumulants for the net-proton distributions are treated as the main observables to search the critical point in experiment.

Recently, the STAR collaboration has measured the energy dependent moments of net-proton multiplicity distributions in Au+Au collisions at = 7.7, 11.5, 19.6, 27, 39, 62.4 and 200 GeV Aggarwal:2010wy (); Adamczyk:2013dal (); Luo:2015ewa (). With the maximum transverse momentum increased from to GeV, the cumulant ratio for net protons shows large deviation from the Poisson expectations and presents an obvious non-monotonic behavior in central Au+Au collisions Luo:2015ewa (), which shows the potential to discover the QCD critical point.

To quantitatively study these experimental measurements, one needs to develop dynamical models near the QCD critical point. Based on linear sigma model, Paech and her collaborators have established the framework of chiral hydrodynamics through coupling the hydrodynamical equations for quarks with the evolution equations for the long wavelength mode of the sigma fieldPaech:2003fe (). Later on, the Frankfurt group further developed chiral hydrodynamics with dissipation and noise Nahrgang:2011mg (); Nahrgang:2011vn (), and with the contributions from the Polyakov loop Herold:2013bi (). Besides the dynamical evolution, an equation of state with a critical point has also been constructed and been applied to pure hydrodynamic simulations Nonaka:2004pg (); Asakawa:2008ti ().

As far as we know, most of the dynamical models near the QCD critical point only
focus on the dynamical evolution of the bulk matter Paech:2003fe (); Nahrgang:2011vn (); Nahrgang:2011mg (); Herold:2013bi (). Without a proper
treatment of the freeze-out procedure, these models can not be used to directly calculate the fluctuations
of produced hadrons as measured in experiment. In this paper, we will introduce
a freeze-out scheme for the dynamical models near the QCD critical point using a modified distribution function.
To avoid the numerical complexities, we will not investigate the dynamical critical phenomena from
a realistic dynamical evolution, but concentrate on studying the static critical fluctuations of emitted classical particles
from a pure hydrodynamic freeze-out surface. We will demonstrate that, with Poisson and Binomial
statistical fluctuation baselines, the static critical fluctuations from our model calculations
can roughly reproduce energy dependent cumulant and of net protons as measured in experiment.
However, the calculated and with both Poisson and Binomial baselines
are always above the experimental data due to the positive contributions
from the static critical fluctuations. In order to qualitatively and quantitatively describe all the related
experimental data, the dynamical critical fluctuations and more realistic non-critical fluctuation baselines should be investigated
in the near future.

## Ii The formalism

.

In this section, we will introduce a freeze-out scheme for the evolution of the bulk matter near the critical point, and then deduce the formulism to calculate the 2nd, 3rd and 4th cumulants for the correlated fluctuations of the classical particles emitted from the hydrodynamic freeze-out surface near with the presence of an external order parameter field. We emphasis that our formalism aims to investigate the static critical phenomena near the QCD critical point, which does not involve the dynamical evolution of the order parameter field. We also neglect the hadronic scatterings and resonance decays below , but leave the investigations of such effects to the future study.

In traditional hydrodynamics, the classical particles emitted from the freezeout surface can be calculated through the Cooper-Frye formula PhysRevD.10.186 (); Kolb:2000sd ():

(1) |

where is the freezeout hyper-surface, is the normal vector on the freezeout surface, and is the distribution function for the classical particles. For Boltzmann statistics, its co-variant form is written as where is the four-momentum of the classical particles (with , is the physics mass of the particles), is the four-velocity of the fluid and is the decoupling temperature on the freeze-out surface.

In the vicinity of the critical point, we assume the classical particles still
satisfy some modified distributions , but with variable
effective masses that fluctuate in position space
through interacting with the correlated fluctuating order parameter
field . For example, the proton interacts with the sigma field through
the coupling,
its effective mass can be written into two parts: . is
the physical mass, is the variable mass induced by the
sigma filed ^{1}^{2}

With the modified distribution function and the variable effective mass, one can, in principle, calculate the critical fluctuations of produced particles through the Cooper-Frye formula in a dynamical model near the critical point, e.g., chiral hydrodynamics. However, such calculations involve event-by-event simulations of the coupled evolution equations for the sigma field and for the bulk matter. To compare with the experimental data, especially, the higher cumulants of net protons, one needs to generate the fluctuating sigma fields that satisfies specific 2-point, 3-point and 4-point correlators. Unfortranately, this has not been numerically realized in current dynamical simulations.

In this article, we will perform event-averaged calculations for the correlated fluctuations of particles emitted
from the hydrodynamic freeze-out surface through expanding the modified equilibrium distribution functions
^{3}

(2) |

where is the traditional equilibrium distribution function: , represents the fluctuation, and is the covariant Lorentz factor. With such expansion, the 2-point, 3-point and 4-point correlators of are expressed as:

(3) | |||||

(4) | |||||

(5) |

note that denotes the expectation values of the connected correlators since disconnected correlators do not contribute to the related cumulants of final produced hadrons.

The above correlators Eqs.(3-5) require detailed expressions for the correlations of the sigma field. To simplify the calculations, we neglect the feedback interactions from particles and only consider the static critical fluctuations of the sigma field. Following Ref. Stephanov:2008qz (); Stephanov:2011pb (), we use the probability distribution with a potential involve cubic and quartic terms to calculate the correlators of the sigma field, which is written as:

(6) |

with

(7) |

With this probability distribution, the correlators of the sigma field at the tree diagram level are:

(8) | |||||

(9) | |||||

(10) | |||||

where the propagator , , is
the mass of the sigma field, and the corresponding correlation length is .

With the above detailed expressions, one can calculate the cumulants of produced hadrons through integrating Eqs.(3-5) on the hydrodynamic freeze-out surface:

(11) | ||||

(12) | ||||

(13) |

In the following calculations, we denote these cumulants of critical fluctuations , and as , and , respectively.

In Eqs.(11-13), the effects from the dynamical evolution of the bulk matter have been partially accounted through these multi-dimensional integrations on the hydrodynamic freeze-out surface. However, our formulism still belongs to the category of static critical fluctuations since the input correlators, Eqs.(8-10), are static correlators without involving a time evolution of the sigma field. If replacing the related integrations of Eqs.(11-13) by the integrations over the whole position space, the standard formula for a static and infinite medium deduced by Stephanov in 2009 Stephanov:2008qz () can be reproduced (please see the appendix for details). Recently, Mukherjee and his collaborators Mukherjee:2015swa () have found that the cumulants of the sigma field could change their sign during the dynamical evolution. However, their method can not be directly implemented to our formulism since the related calculations only consider the zero mode of sigma field, which already erase the needed spacial information.

## Iii Set ups

.

Eqs.(11-13) require a (3+1)-dimensional freeze-out surface for the related multi-dimentional integrations. To obtain such freeze-out surface, we implement the viscous hydrodynamic code VISH2+1 Song:2007fn () and extend its (2+1)-d freeze-out surface to the longitudinal direction within the measured momentum rapidity range through the rapidity correlations between momentum and space Song:2010aq (). The hydrodynamic simulations start at with smooth initial conditions generated from the MC-Glauber model. The normalization factor of the initial entropy density profiles are tuned to fit the multiplicity of pions in central Au+Au collisions at 7.7, 11.5, 19.6, 27, 39, 62.4 and 200 GeV. Since this paper does not aim to fit the flow data at lower collision energies, we directly set the specific shear viscosity and the specific bulk viscosity as once used in Ref. Song:2010mg (); Song:2012ua (). The hydrodynamic freeze-out surface is defined by a constant temperature, which is set to the chemical freeze-out temperature extracted by the statistical model at various collision energies Andronic:2005yp (); Cleymans:2005xv ().

This paper focuses on investigating the fluctuations of net protons. It is thus necessary to nicely fit the related mean values of protons and anti-protons from the hydrodynamic calculations. Since the current version of VISH2+1 still inputs the equation of state s95-PCE with zero chemical potential and does not include the transport equations for the conserved charges, we directly add an effective chemical potential of net baryons to the distribution function in the Cooper-Frye formula, and then fine tune to produce the yields of protons and anti-protons at each collision energies and chosen centralities.

7.7 | 11.5 | 19.6 | 27 | 39 | 62.4 | 200 |

para-I | 3.2 | 2.5 | 2.3 | 2.2 | 2 | 1.8 | 1 | |

6 | 4 | 3 | 2 | 0 | 0 | 0 | ||

14 | 13 | 12 | 11 | 10 | 9 | 8 | ||

1 | 2 | 3 | 3 | 2 | 1 | 0.5 |

para-II | 3.2 | 2.5 | 2.3 | 2.2 | 2 | 1.8 | 1 | |

6 | 4 | 3 | 2 | 2 | 1.5 | 1 | ||

14 | 13 | 12 | 11 | 10 | 9 | 8 | ||

1 | 2.5 | 4 | 4 | 3 | 2 | 1 |

para-III | 2.8 | 1.8 | 1.7 | 1.6 | 1 | 0.5 | 0.1 | |

6 | 4 | 3 | 2 | 2 | 1.5 | 1 | ||

14 | 13 | 12 | 11 | 10 | 9 | 8 | ||

1 | 2 | 3 | 3 | 2 | 1 | 0.5 |

Besides the hydrodynamic freezeout surface, couplings , and the correlation length are additional inputs in our calculations. In order to reproduce the physical mass of protons at zero temperature, in the sigma model Stephanov:2008qz (). At high temperature, calculations from the NJL model show that the coupling decrease as the temperature increase Klevansky:1992qe (). As far as we know, there is no solid calculation on how the coupling changes with temperature and chemical potential along the chemical freeze-out line. We thus treat as a free parameter and tune it within 0 and 10. According to the lattice simulations for the related effective potentials around the critical point, the dimensionless parameters and separately ranges between and from crossover to the first order phase transition Tsypin:1994be (); Stephanov:2008qz (). Considering the critical slowing down near the critical point, the maximum value of the correlation length is set to O(3fm) Berdnikov:1999ph (); Nonaka:2004pg (). Away from the critical point, the correlation length gradually decreases to its natural value around 0.5-1 fm Stephanov:2008qz (); Stephanov:2011pb ().

Table I lists three parameter sets used in our calculations. Here, and are tuned within the above constrains to roughly fit and to describe the decreasing trend of and with Poisson and Binomial baselines. Note that these parameter sets are not unique since the measured cumulants and are always below the Poisson/Binomial baselines, which can not strictly constrain the related parameters of critical fluctuations (please refer to Sec. IV for details).

In order to compare with the experimental data, one also needs to consider
the contributions from trivial statistical fluctuations (noncritical fluctuations).
Here, we assume that the critical and noncritical
fluctuations are independent, which can be treated separately. Following
Ref. Luo:2014tga (), we take either Poisson or Binomial
distributions as the trivial statistical fluctuations. If protons and
anti-protons are independently produced and satisfy the Poisson distributions,
the net protons satisfy the Skellam distributions. The correspondent
cumulants and cumulants ratios are expressed as: , , , , where and are the
mean values of protons and antiprotons. If protons and anti-protons satisfy
independent Binomial distributions, the cumulants of protons and
anti-protons are written as: , ,
,
, where is the mean value of protons/anti-protons and
is an additional parameter, which is determined by . Then, the
cumulants of the net-protons are expressed as: ^{4}

## Iv Numerical results

.

With the formulism and set-ups presented in Sec.II and Sec.III,
we calculate the correlated fluctuations of net protons emitted
from the hydrodynamic freeze-out surface with the presence of an
order parameter field. Fig. 1 and Fig. 2 show the energy dependent
cumulants of net protons in central Au+Au collisions
with either Poisson or Binomial distributions served as the trivial statistical
fluctuation baselines, where the left and right panels
present the results within two different
transverse momentum ranges, and , respectively.
The red stars are the STAR preliminary data Adamczyk:2013dal (); Luo:2015ewa (),
the dashed blue lines are the expectation
values of Poisson/Binomial distributions, and the black and grey curves with
different symbols are the theoretical results from our model calculations.
Here, we implement three parameters sets, as listed in table I, to calculate the
the critical fluctuations ^{5}

In our model calculations, , the mean value of net protons, is used to set the chemical potential of net baryons in the Cooper-fryer formula, which leads to a nice fit of the data within different ranges at various collision energies. are the cumulants to evaluate fluctuations. Generally, the deviations between the data and the Poisson/Binomial baselines are considered as the contributions from the critical fluctuations. After tuning , , and within the allowed parameter ranges, we could roughly describe the decrease trend of and and the non-monotonic behavior of with the increase of collision energy. However, and from our model calculations are always above the Poission/Binomial baselines due to the positive contributions from the static critical fluctuations, which shows even larger deviations from the experimental data. Similarly, the standard critical fluctuations from Stephanov in a static and infinite medium also give positive contributions to and . Such sign problem of and might be solved if further considering the time evolution for the correlated fluctuating order parameter field, which should be studied in the future.

For the energy dependent cumulant , our model calculations can roughly describe the experimental data within both and with the Binomial baselines. However, if using the Poisson baselines, it’s difficult to simultaneously describe within these two different ranges for Au+Au collisions at lower collision energies. At small collision energies below 11.5 GeV, the measured are higher than the poisson expectation values for , but lower than the poisson expectation values for . According to Eqs.(13-15), the change of the ranges only affects the magnitude of of the critical fluctuations, rather than their signs. As a result, our calculations with the Poisson baselines can not simultaneously fit the data within these two ranges.

Fig. 1 and Fig. 2 also show that, with the range changed from to , the measured dramatically increases, showing much larger deviations from the Poison/Binomial baselines. It is thus generally believed that large acceptance is a crucial factor to probe the critical fluctuations in experiment. In our model calculations, the cumulants of the critical fluctuations is closely related to the power of the total net-proton number within the defined and rapidity range. With the maximum increased from to , the averaged number (mean value) of net protons almost increases by a factor of two, leading to a dramatic increase of and in our calculations. One may also notice that, although the correlation length decreased to 1 fm in Au+Au collisions at 7.7 GeV, the critical fluctuations of , and are still large there, especially for the cases with . This is mainly because the net proton numbers within specific acceptance range dramatically increase at lower collision energies. Compared with the case at 19.6 GeV, the mean values of net protons within both ranges almost increase by a factor of 2 for 7.7 GeV, leading to larger , and even with = 1 fm.

Similarly, the higher cumulants of critical fluctuations are dramatically suppressed from central to semi-cental collisions due to the decreased mean values of net protons. To further explore such effects, we extend our calculations to 30-40% centrality bin with the same inputs and parameter sets. Fig. 3 and Fig. 4 show the energy dependent cumulants in 30-40% Au+Au collisions with Poisson and Binomial baselines, respectively. Due to largely reduced critical fluctuations, these black and grey curves from our model calculations are close to the Binomial and Poisson baselines. We have also noticed that the Poisson expectation values are well above the data for different ranges, which make it impossible to fit the if using the Poisson baselines in our model calculations. In contrast, the Binomial baselines are very close to the measured in 30-40% centrality, with which we can roughly fit the data within error bars. However, and are still slightly over-predicted there due to the positive contributions from the critical fluctuations.

Fig. 5 and Fig. 6 show the energy dependent cumulant ratios and in 0-5% and 30-40% Au+Au collisions, where the left and right panels present the results with Poisson and Binomial baselines, respectively. Although and are over-predicted in our model calculations, the cumulant ratios and show better agreement with the experimental data in central Au+Au collisions, except for within (pannel (b) and (f) in Fig. 5). For non-central Au+Au collisions at centrality, the Binomial baselines are very close to the experimental data. Together with the dramatically reduced critical fluctuations, this leads to a nice fit of the experimental data in our model calculations. In contrast, the Poisson baselines are well above the measured for almost all collision energies at the 30-40% centrality bin. Due to small contributions from the critical fluctuations, it is impossible to fit at higher collision energies above 39 GeV in our model calculations with the Poisson baselines.

## V Summary and discussions

.

In this paper, we introduced a freeze-out scheme for the dynamical models near the QCD critical point through coupling the decoupled classical particles with the correlated fluctuating sigma field. With a modified distribution function that satisfies specific static fluctuations, we deduced the formulism to calculate the 2nd, 3rd and 4th cumulants for the particles emitted from the hydrodynamic freeze-out surface with the presence of an order parameter field. We proved that, for a static and infinite medium, this formulism can reproduce the early results of static fluctuations deduced by Stephanov in 2009.

With the parameters tuned within the allowed ranges, we calculated the correlated fluctuations of net protons on the hydrodynamic freeze-out surface in central and non-central Au+Au collisions at various collision energies. With Poisson or Binomial baselines, our model calculations could roughly describe the decrease trend of the cumulants and and the non-monotonic behavior of the cumulant as observed in experiment, but always over-predict the values of and due to the positive contributions from the static critical fluctuations.

After fine tuning the coupling , , and , our model calculations, with the Binomial baselines, can roughly describe the energy dependence of and within different ranges in both central and non central collisions. However, if using the Poisson baselines, it is difficult to simultaneously describe within different ranges in central Au+Au collisions at lower collision energies. Meanwhile, it is impossible to fit and with the much reduced critical fluctuations for Au+Au collisions at 30-40% centrality bin since the Poisson baselines largely deviate from the experimental data there. For our model calculations, the Binomial distributions with two parameters seem to give better non-critical fluctuations baselines when compared with the one-parameter Poisson distributions. However, how the trivial statistical fluctuations, like Poisson distributions, are destroyed through the effects like hadronic scattering and decays, the cut of centrality bins, etc. has not been fully investigated, which should be studied in the near future.

At last, we emphasize that our calculations presented in this paper belong to the category of static critical fluctuations due to the static correlators of the sigma field used in our formulism. In order to qualitatively and quantitatively describe the experimental data, especially to solve the sign problem of and for the critical fluctuations, the dynamical evolution of the order parameter field and other related effects should be further investigated. Along this direction, Mukherjee and his collaborators have done some pioneering work on solving the evolution equations of various cumulants for the sigma field and found that Skewness and Kurtosis could change their sign after the dynamical evolution Mukherjee:2015swa (). However, their approach can not be directly implemented in our formulism on the freeze-out surface since only zero mode of sigma field are considered in their calculations, which already erase the needed spacial information. Dynamical models near the critical point, eg. chiral hydrodynamics, involve a full evolution of the sigma field in position space, together with the coupled evolution of the bulk matter Paech:2003fe (); Nahrgang:2011vn (); Nahrgang:2011mg (); Herold:2013bi (). Combined with the freeze-out scheme introduced in this paper, one could, in principle, calculate the correlated fluctuations of the classical particles decoupled from the bulk matter with the presence of an evolving order parameter field. However, a quantitative study of the higher cumulants of produced hadrons also requires sophisticated simulations for the evolution equation of the sigma field with the properly chosen noise and dissipation terms and initial conditions with correlated fluctuations, etc.. Part of the related work have been done by different groups in the past few years Paech:2003fe (); Nahrgang:2011vn (); Nahrgang:2011mg (); Herold:2013bi (), more progresses are expected in the near future.

## Appendix

.

A. critical fluctuations in static and infinite system.

The standard critical fluctuations in a static and infinite medium were firstly deduced in Stephanov:1999zu (); Stephanov:2008qz () through considering a joint probability distribution for the sigma field and the occupation numbers of the produced particles. In this paper, we assume the coupling between the sigma field and the classical particles introduces a variable effective mass for the particle, which gives a modified distribution function that strongly fluctuates near the critical point. In appendix A, we will show that, for a static and infinite medium, we could re-derive the early results of Stephanov:2008qz () with such modified distribution function.

For a static and infinite medium, integrating Eqs.(3-5) over the whole position space gives:

(14) |

The two point correlator involves an integration of the equal-time propagator of the sigma field:

(15) |

with and , thus

(16) |

where and for a static system.

The 3-point correlator of the sigma field can also be calculated with the variables substitution: , and the Jacobi determination . The integration of the -point correlator can be written as:

(17) | |||||

Then, the 3-point correlator of particle density in momentum space is

(18) |

Following the same steps, we can calculate the 4-particle correlator in momentum space with the variables substitution, which gives

(19) | |||||

With Eqs.(18-19), we now reproduce the static fluctuations for a static and infinite medium
once deduced in paper Stephanov:2008qz ().

B. cumulants of protons anti-protons and net protons.

In Fig. 1-4, we only present the cumulants of net protons for clearness of the figures.
In appendix B, we will compare the cumulants of protons, anti-protons and net protons
in central Au+Au collisions within . Fig. 7
shows the energy dependence of these cumulants with either Poisson or Binomial
baselines. For clearness, we only choose one parameter set (para-I)
to calculate the critical fluctuations. For anti-protons, the cumulants with both critical and non-critical
fluctuations (para-I) are almost overlap with non-critical fluctuation baselines due to the small mean values
at various collision energies which largely suppress the critical fluctuations.
Basically, the cumulants of net protons follows the trends of protons. As already discovered in Fig. 1
and Fig. 2, our calculations can roughly fit for protons and anti-protons with pare-I,
but over-predict and due to the positive contributions from the critical
fluctuations.

C. momentum dependence of critical fluctuations.

In this appendix C, we explore the acceptance dependence of critical fluctuations in the transverse momentum space. The critical fluctuations in our calculations are mainly influenced by two factors: the mean value (average number) of net protons within specific acceptance window and the correlation length . Fig. 8 shows the collision energy dependent cumulants and of net protons, calculated from Eqs.(11-13) within different windows and within momentum rapidity . At larger collision energies above 39 GeV, all cumulants almost approach zero because of the smaller correlation length and reduced mean values of net protons within specific acceptance window. At lower collision energies, increases as the window is broadened. In fact, the mean value of net protons increases with the increase of the -acceptance, leading to the enhanced signals of the critical fluctuations. Note that with the acceptance increased to , almost reach saturations since most of protons and anti-protons are produced below 2 GeV. Fig. 8 also shows that gradually increase as the collision energy decreases. Meanwhile, the correlation length in table I decrease from O(3fm) to 1 fm. As discussed in Sec. V, both correlation length and the mean value of net protons influence the critical fluctuations. At lower collision energy, the drastically increased becomes the dominant factor, which leads to the increasing trend of with the decrease of the collision energy.

In a recent paper Ling:2015yau (), Ling and Stephanov have investigated the acceptance dependence of critical fluctuations in both transverse momentum and rapidity windows, using a simplified freeze-out surface constructed from the blast-wave model. They demonstrated that the window dependence is significant for different order cumulants, which is qualitatively agree with what we find in Fig. 8. They also found that, for typical experimental rapidity acceptance window , larger acceptance leads to significantly increased critical fluctuations. In this appendix, we will not further explore such rapidity acceptance dependence since the 3+1-d freeze-out surfaces used in our calculations are constructed from the 2+1-d freeze-out surfaces of VISH2+1 using the momentum rapidity and space rapidity correlations.

Acknowledgments

We thank M. Asakawa, A. Dumitru, D. Hou, S. Gupta, Y. Liu, X. Luo, M. Stephanov, N. Xu for valuable discussions. This work is supported by the NSFC and the MOST under grant Nos. 11435001 and 2015CB856900.

### Footnotes

- We have shifted the average value of the sigma field to zero. here represents the fluctuations of the sigma field.
- Based on similar idea, Ref. Stephanov:2011pb () has derived the standard critical fluctuations in a static and infinite medium Stephanov:2008qz () with a fluctuating momentum-space distribution function .
- Such equilibrium distribution functions is a stationary solution for the Boltzmann equation
with the presence of an external sigma field Stephanov:2009ra ().
- As pointed in Ref. Kapusta:2011gt (); Young:2014pka (), a hydrodynamic system has imprinted thermal fluctuations according to the fluctuation-dissipation theorem. Considering the numerical complexities, we do not further explore such thermal fluctuations, but only take the Poisson or Binomial distributions as the trivial statistical fluctuations.
- The critical fluctuations are calculated through integrating the static correlators of over the whole freeze-out surface. If constraining the correlations within certain time slice, e.g. , , and at lower collision energies would respectively reduce by O(10%), O(30%) and O(50%) with the same parameter sets as inputs. After re-tuning the related parameters, similar results as shown in Fig. 1 and Fig. 2 could be obtained even with such constraint.

### References

- M. M. Aggarwal et al. [STAR Collaboration], arXiv:1007.2613 [nucl-ex].
- F. Karsch and E. Laermann, In *Hwa, R.C. (ed.) et al.: Quark gluon plasma* 1-59 [hep-lat/0305025].
- Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, Nature 443, 675 (2006).
- Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg and K. K. Szabo, JHEP 0906, 088 (2009).
- A. Bazavov et al. [HotQCD Collaboration], Phys. Rev. D 90, no. 9, 094503 (2014).
- S. P. Klevansky, Rev. Mod. Phys. 64, 649 (1992).
- K. Fukushima, Phys. Lett. B 591, 277 (2004).
- C. D. Roberts and A. G. Williams, Prog. Part. Nucl. Phys. 33, 477 (1994).
- J. Berges, N. Tetradis and C. Wetterich, Phys. Rept. 363, 223 (2002).
- M. A. Stephanov, Prog. Theor. Phys. Suppl. 153, 139 (2004) [Int. J. Mod. Phys. A 20, 4387 (2005)].
- M. A. Stephanov, K. Rajagopal and E. V. Shuryak, Phys. Rev. Lett. 81, 4816 (1998).
- M. A. Stephanov, K. Rajagopal and E. V. Shuryak, Phys. Rev. D 60, 114028 (1999).
- M. A. Stephanov, Phys. Rev. Lett. 102, 032301 (2009).
- C. Athanasiou, K. Rajagopal and M. Stephanov, Phys. Rev. D 82, 074008 (2010).
- Y. Hatta and M. A. Stephanov, Phys. Rev. Lett. 91, 102003 (2003).
- M. Kitazawa and M. Asakawa, Phys. Rev. C 86, 024904 (2012).
- M. M. Aggarwal et al. [STAR Collaboration], Phys. Rev. Lett. 105, 022302 (2010).
- L. Adamczyk et al. [STAR Collaboration], Phys. Rev. Lett. 112, 032302 (2014).
- X. Luo [STAR Collaboration], PoS CPOD 2014, 019 (2014).
- K. Paech, H. Stoecker and A. Dumitru, Phys. Rev. C 68, 044907 (2003).
- M. Nahrgang, S. Leupold, C. Herold and M. Bleicher, Phys. Rev. C 84, 024912 (2011).
- M. Nahrgang, C. Herold, S. Leupold, I. Mishustin and M. Bleicher, J. Phys. G 40, 055108 (2013).
- C. Herold, M. Nahrgang, I. Mishustin and M. Bleicher, Phys. Rev. C 87, no. 1, 014907 (2013).
- C. Nonaka and M. Asakawa, Phys. Rev. C 71, 044904 (2005).
- M. Asakawa, S. A. Bass, B. Muller and C. Nonaka, Phys. Rev. Lett. 101, 122302 (2008).
- F. Cooper and G. Frye, Phys. Rev. D 10, 186 (1974).
- P. F. Kolb, J. Sollfrank and U. W. Heinz, Phys. Rev. C 62, 054909 (2000).
- M. A. Stephanov, Phys. Rev. Lett. 107, 052301 (2011).
- M. A. Stephanov, Phys. Rev. D 81, 054012 (2010).
- S. Mukherjee, R. Venugopalan and Y. Yin, Phys. Rev. C 92, no. 3, 034912 (2015).
- H. Song and U. W. Heinz, Phys. Lett. B 658, 279 (2008); Phys. Rev. C 77, 064901 (2008); H. Song, Ph.D thesis, Ohio State Univeristy (2009), arXiv:0908.3656 [nucl-th] .
- H. Song, S. A. Bass and U. Heinz, Phys. Rev. C 83, 024912 (2011).
- H. Song, S. A. Bass, U. Heinz, T. Hirano and C. Shen, Phys. Rev. Lett. 106, 192301 (2011); Phys. Rev. C 83, 054910 (2011).
- H. Song, Nucl. Phys. A 904-905, 114c (2013), [arXiv:1210.5778 [nucl-th]].
- A. Andronic, P. Braun-Munzinger and J. Stachel, Nucl. Phys. A 772, 167 (2006).
- J. Cleymans, H. Oeschler, K. Redlich and S. Wheaton, Phys. Rev. C 73, 034905 (2006).
- M. M. Tsypin, Phys. Rev. Lett. 73, 2015 (1994); hep-lat/9401034; hep-lat/9601021.
- B. Berdnikov and K. Rajagopal, Phys. Rev. D 61, 105017 (2000).
- J. I. Kapusta, B. Muller and M. Stephanov, Phys. Rev. C 85, 054906 (2012)
- C. Young, J. I. Kapusta, C. Gale, S. Jeon and B. Schenke, Phys. Rev. C 91, no. 4, 044901 (2015)
- X. Luo, B. Mohanty and N. Xu, Nucl. Phys. A 931, 808 (2014).
- B. Ling and M. A. Stephanov, arXiv:1512.09125 [nucl-th].