Non-equilibrium cumulants within model A from crossover to first-order phase transition side

Searching for the location of the critical point and depicting the QCD phase diagram is one of the fundamental topics in heavy-ion physics, which has been extensively studied for decades Klevansky:1992qe ; Roberts:1994dr ; Berges:2000ew ; Karsch:2003jg ; Fukushima:2003fw ; Aoki:2006we ; Fu:2007xc ; Aoki:2009sc ; STAR:2010vob ; Qin:2010nq ; Jiang:2013yoa ; Luo:2017faz . Induced by the divergent fluctuations at the critical point, dramatic increases of the cumulants of the final state proton multiplicities in the critical region are predicted as a consequence of strong coupling between the quark matter and the chiral order parameter field (the $\sigma$ field) Stephanov:1998dy ; Hatta:2003wn ; Stephanov:2008qz ; Athanasiou:2010kw ; Kitazawa:2012at . Even though the correlation length of the chiral order parameter will be suppressed by the finite size of the QCD fireball in heavy-ion collisions, the high-order correlators remain quite sensitive to the increase of the correlation length near the critical point since they are proportional to the high power of the correlation length Stephanov:2008qz . Especially, the quartic critical cumulant is predicted to be negative as approaching the critical point from the crossover side, but positive from the first-order phase transition side, thus the non-monotonic kurtosis along the chemical freeze-out line is expected to carry the signals of the QCD critical point in laboratory observations Stephanov:2008qz ; Asakawa:2009aj ; Friman:2011pf ; Stephanov:2011pb .

The experimental exploration of the QCD phase diagram is performed at the BNL Relativistic Heavy Ion Collider (RHIC). The STAR collaboration has measured the high-order cumulants of net proton production in Au+Au collisions for the collision energy ranging from $7.7$ GeV to $200$ GeV STAR:2010vob ; STAR:2013gus ; Luo:2015ewa . The experimental data of $\kappa\sigma^{2}$ ($=C_{4}/C_{2}$) shows a large deviation from the baselines determined by Poisson statistics and presents complex non-monotonic structure in the collision energy region below $39$ GeV  Luo:2015ewa ; STAR:2020tga , which has not been fully explained in theoretical calculations until now. Meanwhile, more experimental measurements are being executed by the HADES collaboration at GSI HADES:2020wpc and by the ALICE Collaboration at the LHC ALICE:2019nbs . Future experimental facilities such as FAIR in Darmstadt, NICA in Dubna, HIAF in Huizhou and J-PARC in Tokai will further explore the details of the QCD phase transition at the (relative) low collision energy region.

Since the high-order cumulants of the protons are measured after hadrons chemically freeze-out, to identify the signals of the critical point from the experimental data, a proper freeze-out scheme near the critical point is introduced to the hydrodynamic model and the equilibrium cumulants of protons on the freeze-out surface are investigated  Jiang:2015hri ; Jiang:2015cnt ; Jiang:2016duk . The resulted cumulants qualitatively describe the acceptance dependence of the experimental data and roughly fits the $C_{4}$ and $\kappa\sigma^{2}$ data, but $C_{2}$ and $C_{3}$ are overestimated in this framework. Incontrovertibly, the off-equilibrium dynamics play an essential role in the dynamical evolution with phase transition to address the critical slowing down near a critical point Wu:2021xgu ; Berdnikov:1999ph and to reveal the domain formation at the first-order phase transition Randrup:2010ax . Indeed, the high-order cumulants vary significantly, not only in the magnitudes but also the signs, in the real-time evolution compared with the equilibrium hypothesis Mukherjee:2015swa ; Jiang:2017sni ; Jiang:2017fas ; Jiang:2017mji . Furthermore, the early time fluctuations and the critical enhancements around the critical point can be probed by the rapidity-window-dependent Gaussian cumulant  Sakaida:2017rtj . In Ref. Nahrgang:2018afz , it shows that the effects of the nonlinear coupling and finite size are manifested through the reduction of the correlation length near the pseudo-critical temperature. In Ref. Rougemont:2018ivt , the well separated equilibrating time of non-hydrodynamic quasi-normal modes in different channels is investigated at the critical point by a phenomenological realistic holographic model, likely to be model B.

According to the classification of dynamic universality, the QCD chiral phase transition in heavy-ion collisions belongs to Model H, since it is governed by five conserved hydrodynamic modes Rev1977 ; Ma1976 ; Son:2004iv ; Chao:2020kcf ; Nahrgang:2020yxm ; De:2020yyx ; Schweitzer:2021iqk . Dynamical models with hydrodynamic background, such as chiral hydrodynamics Paech:2003fe ; Nahrgang:2011mg ; Nahrgang:2011vn ; Herold:2013bi ; Herold:2014zoa and Hydro+ Stephanov:2017ghc ; Rajagopal:2019xwg ; Du:2020bxp have been developed to describe the medium expansion and the critical mode evolution of QCD phase transition in heavy-ion collisions. However, besides the numerical complexities of these models, the equation of state, the transport coefficients, and the freeze-out scheme in the vicinity of the critical point still needs careful discussions. The numerical calculations of cumulants are even more challenging when the long range correlation of high-order fluctuations An:2020vri and the nonzero-momentum critical modes are considered for different phase transition scenarios, thus it is worth employing the simplified relaxational model A and/or diffusive model B  Mukherjee:2015swa ; Jiang:2017sni ; Jiang:2017fas ; Jiang:2017mji ; Sakaida:2017rtj ; Nahrgang:2018afz ; Rougemont:2018ivt ; Nahrgang:2020yxm to study the dynamical phase transition at a broader chemical potential regions.

Based on the event-by-event Langevin equation, we detailed study the dynamical evolution of the $\sigma$ field’s cumulants in different phase transition scenarios and on the hypothetical freeze-out line. Some preliminary results of the dynamical cumulants in various phase transition scenarios are provided in the proceedings Jiang:2017fas ; Jiang:2017mji . Here, more calculations and discussions are presented, including those involving the impact of the initial temperature on evolution and the dynamical cumulants on the hypothetical freeze-out lines. The paper is organized as follows: In Sec. 2, we introduce the dynamical model in detail where the necessary input parameters are exhibited to solve the stochastic equation, including the initial profile, expansion routine of temperature, damping coefficient, and the magnitude of noise. In Sec. 3, we present the numerical results of $\sigma$ field’s cumulants along the given evolution trajectories in both the crossover and the first-order phase transition scenario. Afterward, we illuminate the results of non-equilibrium cumulants on the hypothetical freeze-out lines. In Sec. 4, we close the article by summarizing our main results and discussing future developments.

Linear sigma model is a QCD-inspired low energy effective theory, depicting the phase structure of strongly interacting matter in the $\mu$-$T$ plane via the order parameter $\sigma=\left\langle\bar{\psi}\psi\right\rangle$ Jungnickel:1995fp ; Skokov:2010sf . As the mass of the $\sigma$ field approaches to zero near the critical point, its correlation length grows to infinity. The corresponding equation of motion of the long wavelength mode is well described by the Langevin equation Nahrgang:2011mg ; Xu:1999aq ; Meistrenko:2020nwx :

where the effective potential of the $\sigma$ field is explicitly written as

$U\left(\sigma\right)$ denotes the vacuum contribution and takes the form of

The parameters $\lambda$, $\sigma$, $h_{q}$, and $U_{0}$ are set by the hadrons properties at zero temperature. As the chiral symmetry is spontaneously broken in the vacuum, the nonzero expectation of the $\sigma$ field is $\left\langle\sigma\right\rangle=f_{\pi}=93$ MeV with $\left\langle\vec{\pi}\right\rangle=0$. In reality, the chiral symmetry is explicitly broken by the light quark mass, then the linear term is included with $h_{q}=f_{\pi}m_{\pi}^{2}$ and $m_{\pi}=138$ MeV. $\nu^{2}=f_{\pi}^{2}-m_{\pi}^{2}/\lambda^{2}$ and the mass of $\sigma$ is $m_{\sigma}\sim 600$ MeV by setting $\lambda^{2}=20$. The zero-point energy $U_{0}=m_{\pi}^{4}/\left(4\lambda^{2}\right)-f_{\pi}^{2}m_{\pi}^{2}$. Note that we have neglected the meson fluctuations of $\vec{\pi}$, since the mass of the triplet is finite in the critical regime. $\Omega_{\bar{q}q}$ represents the contributions from the thermal quarks, which is

$d_{q}=12$ is the degeneracy factor of the quarks. The energy dispersion of the valence quark is $E=\sqrt{p^{2}+m_{q}^{2}(\sigma)}$ with the dynamical quark mass $m_{q}\left(\sigma\right)=m_{0}+g\sigma$ Stephanov:2011pb ; Jiang:2015hri . For $g=3.3$, the mass of the constituent quark is approximately $310$ MeV, and the corresponding proton mass $m_{p}\sim 930$ MeV.

According to the effective potential of Eq. (2), the phase diagram is plotted in Fig. 1 as the function of ($\mu,T$). The dot line denotes crossover at small $\mu$, and the solid line represents the first-order phase transition at large $\mu$. The critical point locates at $(\mu_{cp},T_{cp})\sim\left(205,100.2\right)$ MeV.

In Eq. (1), the damping coefficient $\eta$ and the white noise $\xi\left(t,x\right)$ originate from the interaction between the $\sigma$ field and the heat bath, which is consisted of the thermal quarks Biro:1997va . In this work, $\eta$ is treated as a free parameter, within the range of models allowing, whose values $\eta=1,3,7$ fm^-1 are taken in the following calculations. In the zero-momentum mode limit, the correlation of the noise has the form Nahrgang:2011mg

We remark that the zero-momentum approximation only suits the critical point scenario at the thermodynamic limit. In the realistic case with finite correlation length, we set the spatial noise at different time steps as,

where $G\left(x\right)$ is a random number generator of the standard normal distribution.

For our numerical implementation, we first construct the initial profiles of the $\sigma$ field using the probability function: $P\left[\sigma\right]\sim\exp\left(-E\left(\sigma\right)/T\right)$, with $E\left(\sigma\right)=\int d^{3}x\,\frac{1}{2}\left(\nabla\sigma\left(x\right)\right)^{2}+V_{\mathrm{eff}}\left(\sigma\left(x\right)\right)$. In order to solve Eq. (1), the space-time information of the local temperature, $T(t,x,y,z)$, and baryon chemical potential, $\mu(t,x,y,z)$, have to be known which in principle shall be extracted from the heat bath. For simplicity, we assume the system evolves along the constant baryon density trajectories (seen traj. I and traj. II in Fig. 1), while the spatial-uniform temperature decreases in a Hubble-like way Mukherjee:2015swa :

where $T_{0}$ is the initial temperature, and $t_{0}=1$fm is the initial time. The whole simulation is run in a $V=6.8^{3}$ fm³ box. The space step size $dx=dy=dz=0.2$ fm and the time step size is $dt=0.1$ fm/c. Note that the system volume will affect the configurations of $\sigma$ field in each event, due to the change in the magnitude of the noise term, but has no influence on the results of the event-averaged $\sigma$.

With all the ingredients in hand, we complete the event-by-event simulations of the $\sigma$ field from Eq. (1). In the numerical calculation process, the configurations of the $\sigma$ field at every time step over $10^{5}$ events are recorded. The moments of the $\sigma$ field are then calculated by:

where $\sigma=\frac{1}{V}\int d^{3}x\,\sigma\left(\mathbf{x}\right)$. The non-equilibrium cumulants of the $\sigma$ field are iteratively determined by the following formula:

To summarize, in this work, we study the dynamical evolution of the $\sigma$ field based on the event-by-event simulations of a single component Langevin equation. The temperature decrease for the system is set to be Hubble like and the computation is completed in a finite-size system. We statistically weight the dynamical variable $\sigma$ over $10^{5}$ events during the real-time evolution to obtain its high-order cumulants. With the current model setting, we find that the non-equilibrium cumulants express clear memory effects, and the magnitude of $C_{2}$ slowly increases as the system approaches its critical regions. The signs of $C_{3,4}$, as well as the magnitudes, can differ from the equilibrium ones below $T_{c}$. We also find that the high-order cumulants are significantly enhanced on the boundary side of the first-order phase transition. Finally, the spread of out-of-equilibrium cumulants along the hypothetical freeze-out lines has been presented. The non-equilibrium $C_{3}$ is positive at large baryon chemical potentials, in contrast to the negative sign of the equilibrium $C_{3}$. In the vicinity of CP, with certain parameter sets, the non-equilibrium $C_{4}$ expresses non-monotonic curves at large $\mu$ region. We conclude that the combination of supercooling effect and dynamical effects on the first-order phase transition side plays a dominant role in the nonmonotonicity of the high-order cumulants on the hypothetical freeze-out lines.

Note that the apparent memory effects of the cumulants up to the fourth order based on the Langevin framework are in accord with those obtained by the use of the Fokker-Plank equations Mukherjee:2015swa ; Jiang:2017sni , but currently not sufficiently discussed and presented elsewhere, which is crucial for the study of nonequilibrium signals from experiments. Of course, in order to quantitatively describe the experimentally observed cumulants of net proton multiplicities, more sophisticated and realistic dynamical modeling are required from the theoretical side Kapusta:2011gt ; Akamatsu:2017rdu . A suitable mathematical tool that can accurately describe the properties of fluids in heavy-ion collisions and better performs numerical calculations is still under exploration. It is well known that the method of Fokker-Plank equations is equivalent to the Langevin dynamics in the Markov process Mukherjee:2015swa ; Jiang:2017sni ; An:2021wof . However, taking into account the realistic fireball evolution in heavy-ion collisions, where the spatial distributions of temperature and baryon chemical potential, etc., are nonuniform, the Langevin dynamics has the advantage of easily including those effects in simulations, and further combines with the hydrodynamic equations to describe the complete dynamical process in RHIC. We note here how these non-Markov effects manifest in the Fokker-Plank equations is beyond the scope of this paper. Further investigation and analysis are under work.

In this paper, with a relatively simple setup for the modeling and the parameters, we present the dynamics of the cumulants, which serves as reference information for the analysis of the experimental measurements. For the explanation of the experimental data, a number of effects are playing their roles, including the subject of the proper equation of state Monnai:2019hkn ; Noronha-Hostler:2019ayj ; Parotto:2018pwx , of the unknown parameters of the Ising-to-QCD mapping Mroczek:2020rpm , of the critical transport coefficients Karpenko:2015xea ; Denicol:2015nhu ; Bernhard:2019bmu ; Martinez:2019bsn , of the finite size, finite size scaling and global charge conservation in the vicinity of a CP Sombun:2017bxi ; Antoniou:2017vti ; Vovchenko:2020gne ; Poberezhnyuk:2020ayn ; Tomasik:2021jfd , of the non-critical baselines for the cumulants of net-proton number fluctuations Braun-Munzinger:2020jbk , of the nonuniform temperature/chemical potential effects Zheng:2021pia , and of the proper freezeout scheme in the critical region Pradeep:2022mkf . Besides, further connections between the criticality and the other experimental observable are being established through theoretical efforts Mukherjee:2016kyu ; Sun:2017xrx ; Sun:2018jhg ; Wu:2018twy ; Brewer:2018abr ; Pratt:2019fbj and new technique such as the machine learning method Pang:2016vdc ; Huang:2018fzn ; Jiang:2021gsw is also being developed for the search of QCD phase transition signals as well.

We thank Ulrich. Heinz, Yu-Xin Liu, Swagato Mukherjee, Misha Stephanov, Derek Teaney, Yi Yin, Shanjin Wu and Huichao Song for useful discussion and comments. We thank the anonymous reviewer for his/her valuable suggestions in the impact of spatial non-uniformity of the sigma field. The work of LJ has been supported by the NSFC under grant no. 12105223, and the work of JC has been supported by the start-up funding from Jiangxi Normal University under grant No. 12021211.