Modeling the diffusive dynamics of critical fluctuations near the QCD critical point
MModeling the diffusive dynamics of critical fluctuations near the QCD critical point
Marlene Nahrgang ∗ and Marcus Bluhm † SUBATECH UMR 6457 (IMT Atlantique, Universit´e de Nantes,IN2P3/CNRS), 4 rue Alfred Kastler, 44307 Nantes, France
The experimental search for the QCD critical point by means of relativistic heavy-ion collisionsnecessitates the development of dynamical models of fluctuations. In this work we study the fluc-tuations of the net-baryon density near the critical point. Due to net-baryon number conservationthe correct dynamics is given by the fluid dynamical diffusion equation, which we extend by a whitenoise stochastic term to include intrinsic fluctuations. We quantify finite resolution and finite sizeeffects by comparing our numerical results to analytic expectations for the structure factor and theequal-time correlation function. In small systems the net-baryon number conservation turns out tobe quantitatively and qualitatively important, as it introduces anticorrelations at larger distances.Including nonlinear coupling terms in the form of a Ginzburg-Landau free energy functional weobserve non-Gaussian fluctuations quantified by the excess kurtosis. We study the dynamical prop-erties of the system close to equilibrium, for a sudden quench in temperature and a Hubble-liketemperature evolution. In the real-time dynamical systems we find the important dynamical effectsof critical slowing down, weakening of the extremal value and retardation of the fluctuation signal.In this work we establish a set of general tests, which should be met by any model propagatingfluctuations, including upcoming 3 + 1 dimensional fluctuating fluid dynamics.
I. INTRODUCTION
Conventional fluid dynamics propagates averages ofconserved thermodynamic quantities, like the energydensity or charge densities, requiring approximate lo-cal thermal equilibrium [1]. Small deviations from equi-librium are described by dissipative corrections, whichare quantified by the shear and bulk viscosities and thecharge conductivities or diffusion coefficients. In linearresponse theory these transport coefficients are related tocorrelators of the fluctuations of thermodynamic quanti-ties in the fluid dynamical limit [2, 3]. By the fluctuation-dissipation theorem it is consistent to not only includethe dissipative corrections into the nonlinear fluid dy-namical equations of motion but also the propagation ofthe corresponding intrinsic fluid dynamical fluctuations.These intrinsic fluctuations lead, for example, to non-analytic contributions to the time-dependence of corre-lations [2, 4–8]. But most importantly, they become es-pecially interesting when we study the fluid dynamicalbehavior of a system close to a second-order phase tran-sition [9–11].Developing models and simulations for the real-timedynamics of fluctuations at a phase transition has be-come increasingly important in the field of relativisticheavy-ion collisions. These are performed experimen-tally at the Large Hadron Collider (LHC) at CERN, theRelativistic Heavy-Ion Collider (RHIC) at BNL, the Su-per Proton Synchrotron (SPS) at CERN or the HeavyIon Synchrotron SIS18 at GSI. In the heavy-ion colli-sions strongly interacting matter at extreme tempera-tures T and densities is created [12–14]. The success-ful description of collective effects by conventional fluid ∗ [email protected] † [email protected] dynamical simulations [15–21] and the modification ofhigh-energetic probes measured in heavy-ion collisionscompared to proton-proton collisions [22] are convincingindications for the formation of a new state of matter,the quark-gluon plasma (QGP). At the highest beam en-ergies √ s NN at the LHC the QGP is almost baryon free,i.e. the baryo-chemical potential µ B (cid:39)
0, and the tran-sition to hadronic matter is a crossover as demonstratedby lattice QCD calculations [23]. As the beam energyis lowered, the phase diagram of QCD can be probedat finite net-baryon density [24–28]. An especially in-teresting region in the phase diagram is associated withthe conjectured critical point beyond which the transitionto hadronic matter turns into a first-order phase transi-tion [29–33]. Near the critical point fluctuations in con-served charges are expected to grow large and to imprinton the experimentally observed particle multiplicities inform of large event-by-event fluctuations [34–38]. Indeed,first measurements during the beam energy scan phase Iat RHIC and by the HADES experiment at GSI haveshown interesting features in the kurtosis, a fluctuationmeasure associated with the fourth-order cumulant, ofthe net-proton distribution [39–41]. In thermodynamic,i.e. static and infinite, systems these higher-order cumu-lants are known to be in particular sensitive to the growthof the correlation length of the associated critical fluctu-ations [42–44].Up to this day it is unknown quantitatively how criti-cal fluctuations develop in real-time dynamics. Qualita-tively, dynamical fluctuations of the chiral condensate orthe net-baryon density, as two possible order parameters,have been studied in various works [45–68]. The lack of amore quantitative description is mainly due to the chal-lenges that have to be met when including fluctuations into the standard models of heavy-ion collisions, see [69] fora recent review. For the fluid dynamical description it israther straightforward to include criticality on the level a r X i v : . [ nu c l - t h ] J u l of the equation of state [70–72], but the formulation of al-gorithms to treat intrinsic fluctuations in this frameworkremains a challenge [73–82]. For the microscopic trans-port models, where fluctuations are inherently present,the inclusion of a critical point remains complicated.In this work we study the dynamics of fluctuations ina simpler fluid dynamical model, the diffusion equationin one spatial dimension. Our main intent is to reportthe development of an algorithm, which treats fluctua-tions for the crucial long-wavelength modes reliably, andto present corresponding benchmark tests that should bemet by all future approaches that deal with fluid dynam-ical fluctuations. We focus on the net-baryon density,which in the long-time limit becomes the critical modeassociated with the critical point in QCD. We includethe critical physics in the vicinity of the QCD criticalpoint by a Ginzburg-Landau free energy functional, mo-tivated by the 3D Ising universality class. We then testthe presented algorithm for the linear Gaussian limitsin equilibrium. Here, in particular the static structurefactor and the equal-time correlation function are usefulquantities for probing the dynamics of the fluctuations.We then evaluate the dynamical properties of the system,by looking at the dynamic structure factor in equilibriumfirst. Here, we recover the expected dynamical universal-ity class of model B [9]. Next, we investigate the scenarioof a sudden temperature-quench and finally a Hubble-likeevolution of the temperature. We observe effects of crit-ical slowing down, a weakening and a retardation of themaximal signal. II. DIFFUSIVE DYNAMICS NEAR THE QCDCRITICAL POINT
The equations of relativistic fluid dynamics describethe conservation of energy and momentum and of net-charges via ∂ µ T µν = 0 , (1) ∂ µ N µi = 0 . (2)For our purpose we focus on the non-relativistic evolu-tion of the net-baryon number current N µB = n B u µ + j µB ,where the Navier-Stokes expression for the viscous cur-rent is given by j µB = − Γ T ∆ µν ∂ ν (cid:16) µ B T (cid:17) (3)with ∆ µν = u µ u ν − g µν , fluid velocity u µ and mobilitycoefficient Γ. We consider a system that is decoupledfrom the fluid velocity field which we assume to be space-time independent. In this case we recover the diffusionequation ∂ t n B = Γ T ∇ (cid:16) µ B T (cid:17) (4) for the net-baryon density n B . The diffusive dynamicshappens such as to minimize the free energy in the sys-tem. With the thermodynamic relation µ B = δ F /δn B one obtains the diffusion equation generated by the vari-ation of the free energy functional F for a system of spa-tially homogeneous temperature ∂ t n B = Γ ∇ (cid:18) δ F [ n B ] δn B (cid:19) . (5)Since we are interested in the dynamics of intrinsicfluctuations near the critical point we include a stochasticterm to arrive at the stochastic diffusion equation ∂ t n B = Γ ∇ (cid:18) δ F [ n B ] δn B (cid:19) + (cid:126) ∇ · (cid:126)J , (6)where (cid:126)J is a stochastic current given by (cid:126)J = √ T Γ (cid:126)ζ (7)and (cid:126)ζ is a Gaussian spatio-temporal white noise field withzero mean and unit variance. Fulfilling the fluctuation-dissipation theorem the covariance of the stochastic termguarantees that the long-time equilibrium distribution isgiven by P eq [ n B ] = 1 Z exp (cid:18) −F [ n B ] T (cid:19) , (8)normalized by the partition function Z .We choose the free energy functional near the QCDcritical point to be of the following polynomial form in∆ n B = n B − n c with critical density n c : F [ n B ] = T (cid:90) d x (cid:18) m n c (∆ n B ) + K n c ( ∇ n B ) + λ n c (∆ n B ) + λ n c (∆ n B ) + λ n c (∆ n B ) (cid:19) . (9)We note that the chosen Ginzburg-Landau form forthe critical part of the free energy F may be aug-mented by regular contributions. The coupling coeffi-cients can be calculated through the mapping of the 3-dimensional Ising spin model onto a universal effectivepotential [83, 84]. This determines the dependence ofthese couplings on the thermodynamic correlation length ξ within the given universality class as m = 1 ξ ξ , (10) K = ˜ K/ξ , (11) λ = n c ˜ λ ( ξ/ξ ) − / , (12) λ = n c ˜ λ ( ξ/ξ ) − , (13) λ = n c ˜ λ . (14)In principle, the dimensionless couplings ˜ λ , ˜ λ and ˜ λ have universal values as well, but the uncertainty intranslating the spin variables to the QCD phase diagramleads to rather unknown values for these couplings. Wewill use ˜ λ = 1, ˜ λ = 10 and ˜ λ = 3 in this work. Thisimplies that the temperature dependence of the couplingsis determined entirely by the behavior of ξ establishedthrough a matching to the susceptibility of the Isingmodel scaling equation of state [85]. In the work [83, 84]it turned out to be important to include the λ couplingin order to describe the probability distribution of thefluctuations in the spin model. We therefore include thisterm in our study as well, although in a perturbative ex-pansion in ξ /V with volume V this term is suppressedin the scaling regime [42, 86].As can be seen in Fig. 1, the thermodynamic cor-relation length peaks around T c which we choose as T c = 0 .
15 GeV, while the couplings λ and λ havea minimum at T c . There is a region around T c wherethe nonlinear couplings λ and λ are larger than theGaussian mass parameter m . We expect nonlinear ef-fects to be largest here. The critical net-baryon den-sity n c depends on the location of the critical point andthe equation of state. The net-baryon density at chemi-cal freeze-out as a function of √ s NN was obtained fromstatistical model fits using the Hadron Resonance Gasmodel in [89]. Here, maximal values of n B = 0 . / fm are reached at √ s NN ∼ n B = 5 ρ with ρ = 0 . / fm [90]. In this work, we choose a valueof n c = 1 / (3fm ).The above described setup is in general designed forstudying the diffusion dynamics of critical fluctuationsin three spatial dimensions. The numerical frameworkpresented here focusses on the dynamics restricted toone spatial direction. For this purpose, we scale out thetransverse area A and consider the dynamics only in thelongitudinal direction which resembles the situation metin a highly anisotropic heavy-ion collision. With this, thestochastic diffusion equation Eq. (6) becomes ∂ t n B ( x, t ) = Dn c (cid:0) m ∇ x n B − K ∇ x n B (cid:1) + D ∇ x (cid:18) λ n c (∆ n B ) + λ n c (∆ n B ) + λ n c (∆ n B ) (cid:19) + (cid:112) Dn c /A ∇ x ζ x ( x, t ) , (15)where we have expressed the mobility coefficient Γ = Dn c /T via the diffusion coefficient D and the covariancereads (cid:104) ζ x ( x, t ) , ζ x ( x (cid:48) , t (cid:48) ) (cid:105) = δ ( x − x (cid:48) ) δ ( t − t (cid:48) ). III. EQUILIBRIUM FLUCTUATIONS
In this section we investigate the long-time limit for thestochastic diffusion of the net-baryon density at variousfixed thermal conditions. For this purpose, we consider a c ξ / ξ m fm λ fm λ fm FIG. 1. Scaled temperature dependence of the parametersin the Ginzburg-Landau free energy functional F in Eq. (9).We choose T c = 0 .
15 GeV and ξ = ξ = 0 .
479 fm at T = T = 0 . λ = 1 / fm (notshown) is set constant as a function of T . The temperaturedependence of ξ/ξ , which serves as input for the parametersin this work, follows from a matching to the susceptibility ofthe Ising model scaling equation of state for constant µ B onthe crossover side of the QCD phase diagram, see [87, 88] forsome details. system in a quasi one-dimensional box of length L withperiodic boundary conditions. Initially, the net-baryondensity is constant and set to n B ( x ) = n c . Both, thediscretization with ∆ x = L/N x ( N x is the number ofsites) and the finite size of the box will introduce effectswhich make the results differ from the continuum limit(∆ x →
0) and the thermodynamic limit ( L → ∞ ). Whilethe limited resolution is a technical issue, the finite sizereflects the situation of the fireball created in a heavy-ioncollision. After initialization we let the system equilibrateduring a long time, which is proportional to L /D , beforeevaluating the physical observables such as the varianceand kurtosis or the equal-time corelation function andstructure factor of the system. These are related to theequilibrium distribution which is an invariant measureand independent of D . The latter is exemplarily set to D = 1 fm.We note that the determination of equilibrium results,i.e. the long-time behavior, numerically requires a signif-icant amount of statistics. For dissipation in form of dif-fusion any memory on initial conditions is eventually lostand the fluctuation-dissipation balance guarantees ergod-icity of the system. This implies that ensemble averagescan be either obtained by averaging over multiple sam-ples or equally by averaging over time after performinga sufficient amount of equilibration steps proportional to L / ( D ∆ t ). In this work, the high-statistics equilibriumresults have been obtained by combining both methods.We solve the stochastic diffusion equation Eq. (15) nu-merically within a semi-implicit scheme, where the non-linear terms in ∆ n B are treated explicitly. Charge con-servation is respected with very high precision by impos-ing periodic boundary conditions. More details can befound in Appendix A. A. Static structure factor and equal-timecorrelation function in Gaussian models
The stochastic diffusion equation Eq. (15) contains dif-ferent physics cases. For the Gaussian models the non-linear couplings λ i are equal to zero. In this case, ex-act analytic continuum expressions for prominent phys-ical observables are calculable. One of these representsthe dynamic structure factor S ( k, ω ) for wavevector (cid:126)k and frequency ω . It follows directly from the space-timeFourier transform of the stochastic diffusion equation as S ( k, ω ) ≡ V (cid:104) ∆ˆ n B ( k, ω ) ∆ˆ n ∗ B ( k, ω ) (cid:105) = 2 Dn c k ω + [ Dk ( m + Kk ) /n c ] (16)and entails the dynamical space-time spectrum of thefluctuating net-baryon density. We note that for thespatio-temporal white noise field the dynamic structurefactor is S ζ x ( k, ω ) = L (cid:104) ˆ ζ x ˆ ζ x ∗ (cid:105) = 1. From S ( k, ω ) thespatial spectrum at equal time, i.e. the static structurefactor S ( k ), follows from integration over all ω as S ( k ) = 12 π (cid:90) ∞−∞ S ( k, ω ) dω . (17)The simplest version of a Gaussian model is obtainedwhen ˜ K = ˜ λ = ˜ λ = ˜ λ = 0 in Eq. (9). In this case weare left with the Gaussian mass term which gives rise tothe standard diffusion equation. This model serves as areference and was discussed in detail in [64], where thecorrect numerical implementation of Eq. (15) for this casewas verified. From Eq. (16) the static structure factor for˜ K = 0 follows via Eq. (17) as S ( k ) = n c m , (18)which is independent of the wavevector (cid:126)k . Contraryto simple Euler schemes, the semi-implicit scheme ap-plied in our framework achieves highest accuracy for allwavenumbers independent of the time step ∆ t . As weshow in Appendix B, the corresponding structure factorin discretized space-time S k coincides with Eq. (18) andis therefore independent of the lattice spacing ∆ x . In [64]we verified that this is reproduced in our framework.The version with a term of non-zero ˜ K , which describesa kinetic energy in a Klein-Gordon type action or a sur-face tension in diffusion equations, can still be solved an-alytically in the continuum. In this case, which we will call Gauss+surface model, the static structure factor isgiven by S ( k ) = n c m
11 + Kk /m . (19)Due to the finite surface tension the amplitude of thefluctuations becomes suppressed with increasing k .The numerical results presented in this work have beenobtained for ˜ K = 1 in each of the calculations. For ournumerical framework, the static structure factor for theGauss+surface model in discretized space-time reads (seeAppendix B) S k = n c m
11 + Km ∆ x [1 − cos( k ∆ x )] . (20)With increasing resolution, ∆ x →
0, this result convergesto Eq. (19). In Fig. 2 we show the numerical results forthe static structure factor S k as a function of wavenumber κ for fixed box length L = 10 fm and different resolutionsat two different temperatures. As the considered box isfinite in size and the resolution limited by N x only a finitenumber of modes with discrete κ = kL/ (2 π ) are realized.Our numerical implementation reproduces the analyticexpectations for S k from Eq. (20), thus, resolution effectsare well understood. We note that for a resolution of∆ x = (10 / κ (cid:38) κ (cid:46)
10, which are important for thecritical physics, the continuum limit is reached. Close to T c the amplitude of fluctuations for modes with small κ isincreased compared to temperatures further away while S k is rather independent of T for larger wavenumbers.Another prominent observable is the equal-time cor-relation function of density fluctuations in coordinatespace. In the continuum limit it is defined as the Fouriertransform of the static structure factor S ( k ) in Eq. (17)via (cid:104) ∆ n B ( r )∆ n B (0) (cid:105) = (cid:90) d d k (2 π ) d e i(cid:126)k · (cid:126)r S ( k ) . (21)For the quasi d = 1 dimensional system studied in ourwork the equal-time correlation function of density fluc-tuations in the longitudinal direction is given for theGauss+surface model by (cid:104) ∆ n B ( r )∆ n B (0) (cid:105) = n c Am √ K exp (cid:18) − r m √ K (cid:19) . (22)For ˜ K = 1 we recover the standard relation betweenthe Gaussian mass parameter and the correlation lengthgiven by Klein-Gordon theory. The truly realized corre-lation length in the system depends, however, in generalon the surface tension. The integral of Eq. (22) over dis-tances much larger than the correlation length yields thefull weight of the fluctuation, n c / ( Am ). This is the sameas for the pure Gaussian model with vanishing ˜ K , where S k [f m − ] Gauss+surface modelT = 0.5 GeV N x = 64N x = 128N x = 256N x = 512 0.0001 0.001 0.01 0.1 0 10 20 30 40 50 60 S k [f m − ] κ T = T c N x = 64N x = 128N x = 256N x = 512 FIG. 2. The static structure factor (symbols) as a functionof wavenumber κ in the Gauss+surface model ( λ = λ = λ = 0) for different N x = 64 , , ,
512 and fixed L =10 fm. For both temperatures, T = 0 . T = T c , thetheoretical expectations (solid curves) for the static structurefactor in discretized space-time S k , Eq. (20) with k = 2 πκ/L ,are perfectly reproduced. Because of the Hermitian symmetrybetween κ and N x − κ , S k is symmetric about κ = N x /
2. Withincreasing resolution S k approaches the continuum result S ( k )in Eq. (19). Eq. (22) reduces to (cid:104) ∆ n B ( r )∆ n B (0) (cid:105) = ( n c / ( Am )) δ ( r )and the expected uncorrelated Gaussian limit is recov-ered.In [64], the behavior of (cid:104) ∆ n B ( r )∆ n B (0) (cid:105) for the pureGaussian model was studied numerically. For this modelthe correlation function in discretized space-time is givenby (cid:104) (∆ n B ) j (∆ n B ) l (cid:105) = n c / ( Am ∆ x ) δ jl , where j, l cango over all cells. Accordingly, fluctuations are uncor-related over distances larger than the lattice spacing.In our simulations exact net-baryon number conserva-tion is realized over the entire box of finite length L .This leads to corrections which can analytically be un-derstood by imposing the condition of charge conserva-tion (cid:80) l (cid:104) (∆ n B ) j (∆ n B ) l (cid:105) = 0 for any j , see Appendix C.Correspondingly, the expectation for the equal-time cor- −0.002 0 0.002 0.004 0.006 0.008 0.01 0.012 T=0.5 GeV Gauss+surface model < ( ∆ n B ) j ( ∆ n B ) l > L = 5 fmL = 10 fmL = 20 fmL = 40 fm−0.02 0 0.02 0.04 0.06 0 50 100 150 200 250
T=T c < ( ∆ n B ) j ( ∆ n B ) l > r/ ∆ xL = 5 fmL = 10 fmL = 20 fmL = 40 fm FIG. 3. The equal-time correlation function of density fluctu-ations (symbols) in the Gauss+surface model ( λ = λ = λ = 0) for different L/ fm = 5 , , ,
40, fixed ∆ x =(20 / N x , and a representative A = 1 fm .For both temperatures, T = 0 . T = T c , the nu-merical results are found to perfectly agree with the theo-retical expectations (solid curves) when including the finitesize corrections for exact net-baryon number conservation,cf. Appendix C. The correlation function is symmetric in r = | j − l | ∆ x about r = L/ relation function changes to (cid:104) (∆ n B ) j (∆ n B ) l (cid:105) = n c Am (cid:18) δ jl ∆ x − L (cid:19) , (23)which amounts to a constant negative shift that vanisheswith increasing L for fixed resolution. This behavior wasfound to be perfectly reproduced in the numerics, see [64].For the Gauss+surface model similar considerationscan be made. Numerical results for the equal-time corre-lation function (cid:104) (∆ n B ) j (∆ n B ) l (cid:105) as a function of scaleddistance r/ ∆ x = | j − l | are shown in Fig. 3 for fixedresolution ∆ x and various L at two different tempera-tures. We find that the equal-time correlation functionis shifted to negative values at large distances r demon-strating significant anticorrelations. With increasing boxsize L at fixed resolution the negative shift becomes lesspronounced. This behavior is a consequence of exact net-baryon number conservation, see Appendix C. Taking thelatter into account, cf. Eq. (C4), the corresponding ana-lytic expectations agree well with our numerical results,thus, finite size effects in connection with exact chargeconservation are well under control.For temperatures close to T c , (cid:104) (∆ n B ) j (∆ n B ) l (cid:105) be-comes broader and correlations form over larger distancesas one expects from the continuum expression in Eq. (22).Nonetheless, this depends strongly on the size of thebox and finite-size effects in connection with charge con-servation clearly affect the development of the correla-tions. We note that for the larger systems the equilibra-tion times become very long and increasing computer re-sources are needed to produce equilibrated systems andbuild up the expected long-range correlations. In fact,the tiny deviation between theoretical expectations andnumerical results at large r seen in Fig. 3 at T = T c for L = 40 fm is the result of an insufficient equilibrationbefore evaluating the equal-time correlation function. B. Static structure factor and equal-timecorrelation function in the Ginzburg-Landau model
Let us now study the impact of the nonlinear couplingterms in what we call the Ginzburg-Landau model onthe static structure factor and the equal-time correlationfunction. This is shown in Fig. 4 in comparison with theGauss+surface model results for a system of L = 20 fmwith N x = 256 at T = T c . One observes that the in-fluence of the non-zero λ i is the significant reduction of S k at small wavenumbers κ while it is less important forlarger κ . This reduction of the amplitude of fluctuationsat long wavelengths is also reflected in the developmentof spatial correlations. With non-zero λ i , the equal-timecorrelation function is less broad and long-range correla-tions are suppressed. In addition, correlations at smalldistances are less pronounced which consequently reducesthe quantitative impact of exact charge conservation inthe finite-size system. These effects are found to be lessimportant for T further away from T c .The numerical results of the Gauss+surface model cansuccessfully be described by our analytic expectations indiscretized space-time, see Sec. III A. For the Ginzburg-Landau model, instead, no exact analytic expressions canbe derived to compare the numerics with. We note,however, that the numerical results of the Ginzburg-Landau model on the level of 2-point correlations canformally be described by the analytic expressions of theGauss+surface model but with a modified Gaussian massparameter while K is kept fixed. This effective mass, m eff , is larger than m of the Gauss+surface model forany T . Near T c the relative increase of m eff with respectto m is stronger. For the systems studied in this workwe find no additional ∆ x -dependence in m eff within thestatistical uncertainty. S k [f m − ] κ T=T c Gauss+surface modelGinzburg−Landau model S k [f m − ] κ T=T c Gauss+surface modelGinzburg−Landau model −0.02−0.01 0 0.01 0.02 0.03 0.04 0.05 0.06 0 20 40 60 80 100 120 140
T=T c Gauss+surface modelGinzburg−Landau model < ( ∆ n B ) j ( ∆ n B ) l > r/ ∆ x−0.02−0.01 0 0.01 0.02 0.03 0.04 0.05 0.06 0 20 40 60 80 100 120 140 T=T c Gauss+surface modelGinzburg−Landau model < ( ∆ n B ) j ( ∆ n B ) l > r/ ∆ x FIG. 4. Comparison of the static structure factor (upperpanel) and the equal-time correlation function of density fluc-tuations (lower panel) between Ginzburg-Landau model (cir-cles) and Gauss+surface model (squares) for a system of L = 20 fm, N x = 256 and A = 1 fm at T = T c . Thetheoretical expectations in the Gauss+surface model (solidred curves), see Sec. III A, agree with the numerics. The nu-merical results of the Ginzburg-Landau model can formallybe described by the same analytic expressions when replac-ing m by an effective mass m eff that is fitted to describe S k (dashed blue curves). C. Temperature and system-size dependence of thecorrelation length
The continuum expectation of the equal-time correla-tion function in the Gauss+surface model for an infinitesystem is given in Eq. (22). The numerical results in dis-cretized space resemble this form of an exponential decay.This is also the case when taking non-zero λ i into ac-count. As we have seen in Figs. 3 and 4, net-baryon num-ber conservation in the finite-size system results in a neg-ative offset signalling anticorrelations. Still, an exponen-tial form of the correlation function remains. Therefore,we may fit the numerical results of the Gauss+surface ξ ~ / ξ L [fm]
T=T c T=0.5 GeVGauss+surface modelGinzburg−Landau model ξ ~ / ξ L [fm]
T=T c T=0.5 GeVGauss+surface modelGinzburg−Landau model
FIG. 5. L -dependence of the correlation length ˜ ξ in unitsof ξ in the Gauss+surface (squares) and Ginzburg-Landau(circles) models at fixed resolution ∆ x = (20 / T = 0 . T = T c (solid symbols). The horizontal, grey dotted linesshow for comparison the corresponding scaled thermodynamiccorrelation length ξ/ξ for an infinite system, cf. the inputparameters in Fig. 1. and Ginzburg-Landau models with an ansatz that con-tains the exponential behaviour and the offset (see Ap-pendix D for details) in order to determine the correlationlength ˜ ξ . The latter depends besides T in particular onthe system size L and can be different from the thermo-dynamic correlation length ξ .In Fig. 5 we show the system-size dependence of the fit-ted ˜ ξ in the Gauss+surface and Ginzburg-Landau modelsfor two different T at fixed resolution ∆ x = (20 / x -dependence can be estimated to be onthe per cent level for all T and L . For the parametersstudied in this work, cf. Fig. 1, the maximally reachedthermodynamic correlation length in an infinite system, ξ , is about 3 fm near T c and minimally we have ξ = ξ at T = 0 . T = 0 . L = 5 fm is already sufficientfor ˜ ξ to reach approximately the value of ξ . This remainsunchanged with increasing L . However, for all other T with a larger ξ charge conservation turns out to be im-portant, particularly in the smaller systems. In fact, itcan lead to a sizeable reduction of ˜ ξ compared to ξ for L = 5 fm. This effect is pronounced strongest at T c (solidsquares and circles). For T = T c the fitted correlationlength increases strongly toward ξ with increasing L forthe Gauss+surface model. For L = 40 fm one finds ˜ ξ tobe approximately ξ . In contrast, in the Ginzburg-Landaumodel ˜ ξ remains always small compared to ξ and showswithin the statistics a negligible system-size dependencefor L ≥
10 fm. This reduction is entirely a consequenceof the nonlinear interactions. ξ ~ / ξ T/T c ξ / ξ Gauss+surface modelGinzburg−Landau model ξ ~ / ξ T/T c ξ / ξ Gauss+surface modelGinzburg−Landau model
FIG. 6. Comparison of the scaled temperature dependence ofthe correlation length ˜ ξ in units of ξ in the Gauss+surface(squares) and Ginzburg-Landau (circles) models with thescaled thermodynamic correlation length ξ/ξ for an infi-nite system (solid line) for simulations with L = 20 fm and∆ x = L/ In Fig. 6 we compare for L = 20 fm the fitted corre-lation length as a function of temperature with ξ . Oneobserves that ˜ ξ is approximately ξ in the Gauss+surfacemodel for all T except very close to T c , where finite-sizeand charge-conservation effects are strongest, cf. Fig. 5.From this observation we conclude that in order to drawphysical conclusions a reasonable compromise betweenfinite-resolution and finite-size effects on the one handand limited computational resources on the other handis to study systems of L = 20 fm and N x = 256 inthis work. The presence of the nonlinear interactionsin the Ginzburg-Landau model impacts the developmentof long-range correlations significantly. For all T we finda ˜ ξ which is smaller in the Ginzburg-Landau model thanin the Gauss+surface model. While far away from T c the effect is tiny, the reduction is visible in the vicinityof T c . This behaviour is in line with the temperaturedependence of the parameters, see Fig. 1, and with theobservation that for describing the structure factor andthe correlation function in the Ginzburg-Landau modelby the analytic expressions of the Gauss+surface modelone needs m eff > m . In fact, we find that m eff /m be-haves approximately like the ratio of the fitted correlationlengths in the Gauss+surface to the Ginzburg-Landaumodel. We expect that the fluctuation observables aresimilarly affected by this. D. Temperature and system-size dependence ofGaussian and non-Gaussian fluctuations
We now turn to the study of fluctuation observables inthe Gauss+surface and Ginzburg-Landau models. Wewill concentrate on the discussion of local quantities,i.e. on the fluctuations in the net-baryon density con-tained within one grid spacing, on an event-by-event ba-sis. The local variance, σ , is equivalent to the equal-timecorrelation function (cid:104) (∆ n B ) (cid:105) at r = 0. From Eq. (22)we see that σ ∼ ξ . Since the Gaussian mass parameter m ∼ /ξ drops rapidly around T c with a minimum at T c ,cf. Fig. 1, we expect that the local variance is largest at T c in both the Gauss+surface and the Ginzburg-Landaumodel. The local excess kurtosis, κ , is defined as κ = µ σ − , (24)where µ = (cid:104) (∆ n B ) (cid:105) at r = 0 is the fourth central mo-ment of local fluctuations. The excess kurtosis must van-ish for the Gaussian models while in the presence of non-linear coupling terms it provides a measure for the non-Gaussianity of the equilibrium distribution. The localskewness was found to be subject to large statistical un-certainties in the studied finite-size systems with chargeconservation and as a consequence will not be discussedin this work.In Fig. 7 we show numerical results for the system-size dependence of σ and κ in the Gauss+surface andGinzburg-Landau models for two different T at fixed res-olution ∆ x = (20 / σ for an infinite system isgiven by σ = n c A √ Km , (25)which is indicated by the grey dotted lines. In a finite-size system the local variance can be significantly smallerdue to charge conservation, cf. Fig. 3, but increases withincreasing L approaching the limit Eq. (25). The ob-served reduction of σ in the Ginzburg-Landau model isin line with the behavior seen in m eff and ˜ ξ , see Fig. 5and the discussion in section III B. We find a negligibleresidual ∆ x -dependence in σ for all T and L similar to ˜ ξ .This is in contrast to the behavior noted in [64] for thepure Gaussian model where the local variance dependsexplicitly on the resolution, cf. Eq. (23). This unphysi-cal behavior is cured by the inclusion of a finite surfacetension, see also the discussion in [68]. The local excesskurtosis vanishes within the statistical uncertainty in theGauss+surface model. In the Ginzburg-Landau model,instead, κ is non-zero and found to increase in magnitudewith L but also seems to approach a limiting value withincreasing system-size. The residual ∆ x -dependence is abit stronger than for σ but still on the few-percent level.Both σ and κ are significantly larger at T = T c (solidsquares and circles) than at T = 0 . T c finite-sizeeffects in both observables are clearly more pronouncedthan at T = 0 . κ . σ Gauss+surface modelGinzburg−Landau modelT=T c T=0.5 GeV σ Gauss+surface modelGinzburg−Landau modelT=T c T=0.5 GeV −0.4−0.3−0.2−0.1 0 5 10 20 40 κ L [fm]
T=T c T=0.5 GeV
FIG. 7. Results for the local variance σ and local excesskurtosis κ in the Gauss+surface (squares) and Ginzburg-Landau (circles) models for different system sizes L at fixedresolution ∆ x = (20 / T = 0 . T = T c (solid symbols). Thehorizontal, grey dotted lines in the upper panel show for com-parison the corresponding continuum expectations for σ inthe Gauss+surface model for an infinite system, cf. Eq. (25). In Fig. 8 we compare the temperature dependenceof the local variance and local excess kurtosis in theGauss+surface (squares) and Ginzburg-Landau (circles)models for L = 20 fm with N x = 256. The reduction seenin σ for the Ginzburg-Landau model compared to theGauss+surface model is in line with the findings for thetemperature dependence of the fitted correlation lengthin Fig. 6. In fact, within the numerics we find that σ scales approximately as σ ∼ ˜ ξ for all T as expected fromEq. (25). The numerical results for the local excess kur-tosis highlight an important difference between the twomodels: while κ vanishes within the acquired statisticsin the Gauss+surface model, it is non-zero and nega-tive for the chosen values of ˜ λ i in the Ginzburg-Landaumodel. One observes a non-monotonic temperature de-pendence with a prominent peak structure in the vicinityof T c , where ˜ λ and ˜ λ become the dominant parameters,cf. Fig. 1. −0.3−0.2−0.1 0 0.5 1 1.5 2 2.5 3 3.5 κ T/T c σ Gauss+surface modelGinzburg−Landau model
FIG. 8. Results for the local variance σ and local excesskurtosis κ as functions of the scaled temperature T /T c in theGauss+surface (squares) and Ginzburg-Landau (circles) mod-els for simulations with L = 20 fm and ∆ x = L/ IV. DYNAMICS OF GAUSSIAN ANDNON-GAUSSIAN FLUCTUATIONS
We now turn to the study of the dynamics of the sys-tem, which we discuss in three steps: first, we investigatethe dynamical properties in equilibrium in form of thedynamic structure factor, next we study the response ofthe system to a sudden quench in temperature and finallylook at a Hubble-like reduction of the temperature as afunction of time. Note that the dynamical properties de-pend on the value and/or the temporal behavior of thediffusion coefficient D , which as a function of tempera-ture is defined as D = Γ T /n c , where we fix D ( T ) = 1 fmat T = T = 0 . A. Dynamic structure factor and relaxation time
The dynamical properties of the system in equilibriumat a fixed temperature are encoded in the dynamic struc-ture factor. The time-dependence of the spatial spec-trum of the fluctuating net-baryon density is related tothe spatial Fourier transform of the stochastic diffusion S k ,t t [fm] T=T c Gauss+surface modelGinzburg−Landau model S k ,t t [fm] T=T c Gauss+surface modelGinzburg−Landau model
FIG. 9. The dynamic structure factor S k,t as a function oftime for κ = 2 in the Gauss+surface and Ginzburg-Landaumodels for L = 20 fm and N x = 256 at T = T c . equation, Eq. (15), and can be obtained from S ( k, ω ) bythe Fourier transformation into the time-domain viz S ( k, t ) ≡ V (cid:104) ∆ˆ n B ( k, t (cid:48) ) ∆ˆ n ∗ B ( k, t (cid:48) + t ) (cid:105) = 12 π (cid:90) ∞−∞ S ( k, ω ) e iωt dω . (26)For the Gaussian models with S ( k, ω ) given in Eq. (16)this amounts to S ( k, t ) = S ( k ) exp ( − t/τ k ) (27)in the continuum limit, where the static structure fac-tor S ( k ) is given by Eqs. (18) or (19) and the inverserelaxation time reads τ − k = Dm n c (cid:18) Km k (cid:19) k . (28)By setting K = 0 we find the expression of τ k for thepure Gaussian model.Numerically, we study the dynamic structure factor S k,t in discretized space-time by analyzing the correla-tor of the density fluctuations in the mixed representa-tion for modes with given wavevector (cid:126)k and wavenumber κ = kL/ (2 π ), see Appendix E. Exemplarily for κ = 2,we contrast S k,t at T = T c for the Gauss+surface andGinzburg-Landau models in Fig. 9. One clearly observesan exponential decay of the correlator in both modelssimilar to the expected behavior in the continuum limit.As for the static observables, the nonlinear interactions inthe Ginzburg-Landau model reduce the dynamic struc-ture factor compared to the Gauss+surface model and,in addition, accelerate its exponential decay. We notethat in the pure Gaussian model S k,t for the same κ ismuch larger and relaxes significantly slower than in theother models.0 τ k [f m ] T=0.5 GeVGauss+surface modelGinzburg−Landau model 0.001 0.01 0.1 1 τ k [f m ] T=0.5 GeVGauss+surface modelGinzburg−Landau model 0.001 0.01 0.1 1 10 100 0 2 4 6 8 10 τ k [f m ] κ T=T c Gauss+surface modelGinzburg−Landau model 0.001 0.01 0.1 1 10 100 0 2 4 6 8 10 τ k [f m ] κ T=T c Gauss+surface modelGinzburg−Landau model
FIG. 10. Relaxation time τ k (symbols) as a function of κ for L = 20 fm and N x = 256 at T = 0 . T = T c (lower panel) in the Gauss+surface and Ginzburg-Landau models. The results are compared with the contin-uum expression of the Gauss+surface model, Eq. (28), shownas red solid lines, and a modified expression with m replacedby m eff , see Fig. 4, shown as blue dashed lines. The relaxation time in the Gauss+surface model for aspecific mode k can be determined by fitting the corre-sponding S k,t with an exponential ansatz of the form ofthe continuum expression. For T = 0 . T = T c the results for not too large κ are shown in Fig. 10 (redsquares). As one would expect, τ k is drastically enhancednear T c and long-wavelength (small κ ) modes relax sig-nificantly slower than short-wave (large κ ) fluctuations.The continuum results based on Eq. (28) are also shownas red solid lines in Fig. 10. We find that the results of thefits to the data from simulations with ∆ x = (20 / κ (see the discussion in Appendix E).The exponential decay of S k,t seen in Fig. 9 for theGinzburg-Landau model suggests to use a similar ansatzto determine τ k in this case. The results are shown byblue circles in Fig. 10. The nonlinear interactions arefound to reduce the fitted relaxation time, in particular,for modes with small κ , and the effect is more promi-nent in the vicinity of T c . For larger κ , fluctuations areless affected by the nonlinear interactions and τ k in theGauss+surface and the Ginzburg-Landau model is com- parable. The k -dependence of our numerical results for τ k in the Ginzburg-Landau model can quite accuratelybe described by the continuum expression Eq. (28) ofthe Gauss+surface model by replacing m with m eff , seeblue dashed lines in Fig. 10. The values for the modifiedGaussian mass parameter m eff > m are those necessaryfor describing the behavior of the static structure factorin the Ginzburg-Landau model discussed in section III B.The comparison of the fit results with the analytic ex-pectations in the Gauss+surface model indicates that thesimulations carried out with N x = 256 at L = 20 fm arealready sufficiently close to the continuum limit, also forthe dynamic observables. To test further how well an-alytic expectations for resolution effects are reproducednumerically, we decrease the resolution in the simulationsby a factor 4 and consider N x = 64. From Eqs. (E9) -(E11) one expects that a decrease in resolution resultsin an increase of the fitted relaxation time, in particu-lar for larger κ . This is precisely observed in the resultsdepicted in Fig. 11. In fact, the comparison of the fit re-sults for τ k with the expectations for the relaxation time,Eqs. (E9) - (E11), shows that resolution effects are wellcontrolled.The determination of the dynamic structure factor andof the relaxation time allows us to study the correla-tion length dependence of τ k for modes which are cor-related over the distance ˜ ξ . For this purpose, we ana-lyze τ ∗ , the relaxation time of modes with k ∗ = 1 / ˜ ξ ( T ),as a function of T , where ˜ ξ ( T ) is the fitted correla-tion length discussed in section III C. Results for theGauss+surface and Ginzburg-Landau models are shownin Fig. 12 (symbols). We find τ ∗ to behave like a ˜ ξ z withproportionality factor a and dynamic scaling (critical)exponent z . For both models, the best fit (filled bandsin Fig. 12) gives z = 4 ± . a (cid:39) . / ( D fm z − ).This proportionality factor confirms our expectations, a = n c ξ / ( D (1 + ˜ K )), based on the continuum expres-sion of τ k in the Gauss+surface model. We also indicatethat other scaling exponents, e.g. z = 3 (dashed lines)or z = 5 (dotted lines), fail to describe the numericallyrealized scaling with the correlation length. This showsthat our simulations reproduce the dynamic scaling be-havior one would expect for the stochastic diffusion of aconserved charge which is the one of model B within theclassification scheme [9]. B. Relaxation of fluctuation observables after atemperature-quench
The relaxation dynamics of fluctuation observablessuch as the local variance σ and the local excess kur-tosis κ toward equilibrium can be studied through thesudden quench in temperature from a well prepared ini-tial condition. For this purpose, we first let the systemequilibrate at T = T = 0 . τ = τ to threedistinct values T ∗ , namely T ∗ = T c , T ∗ = 0 .
18 GeV1 τ k [f m ] T=0.5 GeVN x = 64N x = 256 τ k [f m ] T=0.5 GeVN x = 64N x = 256 τ k [f m ] κ T=T c N x = 64N x = 256Gauss+surface model τ k [f m ] κ T=T c N x = 64N x = 256Gauss+surface model FIG. 11. Dependence of the relaxation time (symbols) on theresolution in the Gauss+surface model ( λ = λ = λ = 0).We contrast for T = 0 . T = T c (lowerpanel) simulations for N x = 64 (squares) and N x = 256 (tri-angles) at fixed L = 20 fm. Analytic expectations based onEqs. (E9) - (E11) are shown as dashed and solid lines, respec-tively. We note that the analytic results based on Eq. (28)and on Eqs. (E9) - (E11) are practically indistinguishable for∆ x = (20 / κ . and T ∗ = 0 . σ in theGauss+surface model.The results for the relaxation behavior of σ and κ are shown in Fig. 13. One observes that the relaxationdynamics is quite abrupt initially. We find that withdecreasing quench-temperature T ∗ the time it takes σ and κ to relax to the corresponding equilibrium result(horizontal, grey dotted lines) increases. This is to beexpected since for smaller T ∗ we have a smaller diffusioncoefficient D and, moreover, the difference between theequilibrium values at T and at a T ∗ close to T c is larger.In addition, higher moments appear to approach theirequilibrium expectations slower. By increasing the initialvalue of the diffusion coefficient to D ( T = T ) = 2 fm,the relaxation rate is overall increased and the fluctuationobservables relax quicker toward equilibrium, see also thediscussions in [64, 68]. We note that the determinationof the relaxation time of fluctuation observables within a τ * [f m ] T/T c Gauss+surface modelGinzburg−Landau model τ * [f m ] T/T c Gauss+surface modelGinzburg−Landau model
FIG. 12. Scaling behavior of the relaxation time τ ∗ for modeswith k ∗ = 1 / ˜ ξ with the correlation length ˜ ξ in Gauss+surfaceand Ginzburg-Landau models as a function of T /T c . Fitresults to simulations with L = 20 fm and N x = 256(symbols) are contrasted with the scaling ansatz a ˜ ξ z , where a (cid:39) . / ( D fm z − ) and z = 4 ± . z = 3 (dashed lines) and z = 5 (dotted lines) arealso indicated. σ T* = 0.18 GeVT* = T c T* = 0.2 GeV−0.3−0.2−0.1 0 0 2 4 6 8 10 12 14 16 18 κ τ − τ [fm]Ginzburg−Landau model FIG. 13. Relaxation dynamics of the local variance σ andthe local excess kurtosis κ toward their equilibrium valuesunder a sudden quench in temperature from equilibrium at T = 0 . T ∗ at time τ = τ . The shown results are forthe Ginzburg-Landau model with L = 20 fm and N x = 256. C. Time-evolution of fluctuations in a coolingsystem
Assuming a dynamical evolution of the temperatureof the system allows us to highlight some importantnonequilibrium effects. For this purpose, we consider asimple, spatially homogeneous time-dependence of T inthe Hubble-like form T ( τ ) = T (cid:16) τ τ (cid:17) dc s (29)with dimension d = 3, speed of sound c s = 1 / T = 0 . τ = 1 fm at which thesystem is in equilibrium. For this cooling scenario thecritical temperature is reached at τ c − τ = 2 .
33 fm.The time-dependent temperature translates into time-dependent couplings via Eqs. (10) - (14), which are shownin Fig. 14. Due to the fast initial decrease of T and theslower decrease at later times in Eq. (29) the thermody-namic correlation length is more symmetric between theearly and late times than it is in comparison to the highand low temperatures in Fig. 1. Still, all couplings except λ , which is independent of the correlation length, have adip at the time when the critical temperature is reached.This is the region where we expect nonequilibrium effectsto be most prominent.We first study the impact of the dynamical evolutionof T on the equal-time correlation function and the as-sociated correlation length ˜ ξ . The results presented hereare obtained for the Ginzburg-Landau model where weanalyzed a sufficient amount of events. The form of theequal-time correlation function is clearly affected by thedynamics, see upper panel in Fig. 15 (squares) for T = T c .On the quantitative level, this is also determined by thetemporal evolution of the diffusion coefficient D . For nottoo large initial values (such as D ( T ) = 1 fm) it alreadysignificantly decreased (to D ( T c ) = 0 . T c is reached and, thus, local fluctuationscannot rapidly enough be balanced throughout the en-tire finite-size system. As a consequence, correlations atzero distance do not build up quickly enough from thesmaller value at T toward the equilibrium value at T c (see upper panel in Fig. 15 (circles) and Fig. 8) and lagbehind. Around these local fluctuations, anticorrelationsare present due to charge conservation. In the dynamicalsituation they do not have sufficient time to diffuse intothe entire system. We therefore see a dip of the correla-tion function around r = 40 ∆ x , while it approaches zeroat larger distances. This local balancing of the fluctua-tions reduces the correlation length, as we discuss in the τ − τ [fm] ξ / ξ m fm λ fm λ fm FIG. 14. Time-dependence of the parameters entering theGinzburg-Landau free energy functional F in Eq. (9), seeFig. 1 for comparison. The coupling λ = 1 / fm (not shown)is set constant. The critical temperature T c = 0 .
15 GeV isreached at τ c − τ = 2 .
33 fm. following.While the form of the equal-time correlation functionin the dynamical scenario is not in one-to-one correspon-dence with the equilibrium form, we may still analyzethe visible exponential decrease in the region from smallto intermediate distances to deduce a correlation length.By using the ansatz employed in section III C for differ-ent T , i.e. at different times τ − τ , we obtain the resultfor the dynamical ˜ ξ shown in the lower panel of Fig. 15(squares). In comparison to the equilibrium result (cir-cles, cf. also Fig. 6) we observe clear deviations highlight-ing two distinct nonequilibrium effects: first, the overallmagnitude of ˜ ξ is significantly reduced as a consequenceof the dynamics. Secondly, there are clear indications fora retardation effect due to the rapid cooling in T . Thedynamical ˜ ξ remains initially smaller than its equilibriumcounterpart for given T but then develops a maximum ata temperature far below T c such that at late times it isactually larger than in the equilibrium situation. Thepronounced structure traditionally associated with thephase transition is shifted to later times and, thus, dif-ferent thermal conditions. We expect similar effects forthe fluctuation observables.In Fig. 16 we show the temporal evolution of the localvariance σ (upper panel) and the local excess kurtosis κ (lower panel) as a function of time in the Gauss+surface(red bands) and Ginzburg-Landau (blue bands) models.For both models, σ peaks at a time τ − τ shortly af-ter T c is reached during the evolution. The retardationshift appears slightly larger in the Gauss+surface thanin the Ginzburg-Landau model. As in the equilibriumsituation, σ in the Ginzburg-Landau model stays be-low the Gauss+surface model result, but the reduction3 −0.0100.010.020.03 0 20 40 60 80 100 120 T=T c Ginzburg−Landau modeldynamicsequilibrium < ( ∆ n B ) j ( ∆ n B ) l > r/ ∆ x−0.0100.010.020.03 0 20 40 60 80 100 120 T=T c Ginzburg−Landau modeldynamicsequilibrium < ( ∆ n B ) j ( ∆ n B ) l > r/ ∆ x ξ ~ [f m ] τ − τ [fm] FIG. 15. Dynamical behavior (squares) of the equal-timecorrelation function (upper panel) and the fitted correlationlength (lower panel) in a cooling system, cf. Eq. (29), in com-parison with the equilibrium situation (circles) at the same T . The results are obtained for the Ginzburg-Landau modelwith L = 20 fm and N x = 256. We show a snapshot of thecorrelation function as a function of distance r = | j − l | ∆ x with ∆ x = L/N x at time τ c − τ = 2 .
33 fm, i.e. at T = T c ,and ˜ ξ as a function of τ − τ , where the vertical line indicates τ c − τ . of its maximal value due to the dynamics is significantlystronger in the Gauss+surface model (by 46% comparedto 17% for the Ginzburg-Landau model).For the local excess kurtosis we note that in the ab-sence of nonlinear coupling terms κ vanishes in the dy-namical scenario as it did in equilibrium (see red bandand open squares in the lower panel of Fig. 16). For theGinzburg-Landau model κ starts at its equilibrium valuefor T and initially follows the equilibrium behavior forgiven T (see blue band and open circles in the lower panelof Fig. 16). However, within the band of statistical un-certainties it quickly lags behind the equilibrium situa-tion as reflected in the reduced magnitude of κ . We canclearly see that in the dynamical scenario the minimum σ Ginzburg−Landau modelGauss+surface model−0.3−0.2−0.1 0 0 1 2 3 4 5 6 7 8 κ τ − τ [fm] FIG. 16. Dynamical behavior (colored bands) of the lo-cal variance σ and the local excess kurtosis κ for theGauss+surface and Ginzburg-Landau models in a cooling sys-tem, cf. Eq. (29), as a function of time. For the local ex-cess kurtosis we compare with the equilibrium results (opensquares and circles) at the corresponding T . The shown re-sults are for L = 20 fm and N x = 256. in the local excess kurtosis is shifted to a later time than τ c − τ and that the magnitude of this minimum is sig-nificantly reduced (by approximately 30%) compared tothe equilibrium result. At later times, the retardationeffect leads to a dynamical κ slightly larger in magni-tude than in equilibrium. As is evident from Fig. 16, thenonequilibrium effects influence κ stronger than σ in theGinzburg-Landau model. The behavior seen in the fluc-tuation observables resembles qualitatively the one dis-cussed for ˜ ξ above with an important difference: the shiftof the maximum in ˜ ξ to smaller T is larger than in σ or κ . In future work we will investigate the fluctuationobservables over larger subvolumes to see if the relationwith the correlation length is restored. We note that anoverall reduction of the dynamical diffusion coefficient(by lowering its initial value D ( T )) results in a strongerretardation and a stronger reduction of the magnitude ofthe fluctuation signal.4 V. DISCUSSION
In this work we presented a first rigorous implemen-tation of the one-dimensional stochastic diffusion equa-tion near the QCD critical point. First, we benchmarkedthis implementation in the linear approximation, includ-ing Gaussian mass and surface tension terms, versus an-alytic results of the equal-time correlation function andthe static structure factor. Based on these tests, we chosethe resolution of the spatial discretization as to reproducethe behavior of the first 50% of the wavenumbers in thecontinuum limit. Charge conservation is found to playan important role for the correlation function and limitsthe growth of long-range correlations. In the same sensethe growth of the correlation length near T c is limitedfor system sizes up to a few times the thermodynamiccorrelation length.In equilibrium we investigated the temperature depen-dence of the local variance and excess kurtosis. The lattertakes non-zero values as soon as the nonlinear couplingterms in the Ginzburg-Landau free energy functional areincluded. The expected non-monotonic behavior around T c is clearly observed. The inclusion of the nonlinearcoupling terms reduces the variance of the system by afactor of two near T c .From the dynamic structure factor we obtained therelaxation time of the critical mode. It is found toscale with the correlation length according to model Bof dynamical universality. Finally we investigated the re-sponse of the system to changes of the temperature, first,via a sudden quench in temperature and, second, via aHubble-like time evolution. We observe again that thegrowth of the correlation length is limited by charge con-servation effects, this time in a dynamical setup. Here,fluctuations do not have enough time to diffuse to largerdistances and, thus, the correlation length is limited toa smaller range. Fluctuation observables are reduced inmagnitude and shifted to smaller temperatures due tononequilibrium effects. Higher-order cumulants are im-pacted stronger than the variance by the nonequilibriumsituation, i.e. they need more time to relax, the magni-tude of their extrema is more reduced compared to theequilibrium values and the retardation effect is stronger.We emphasize in particular the importance of bench-marking the approach to the dynamics of fluctuationsagainst analytic results, like the correlation function, thestatic and the dynamic structure factor. This shouldbe a standard requirement for all models dealing withthe dynamics of fluctuations, including more complex ap-proaches to fluctuating fluid dynamics.The presented cumulants are evaluated as local ob-servables over individual cells of the simulation region,which serves well the purpose of understanding the basicdynamics of fluctuations in stochastic partial differentialequations. For aiming at a comparison with experimen-tal data from heavy-ion collisions, integrated observablesin finite kinematic regions are of additional interest. Astudy of fluctuations over larger subregions of observa- tion similar to [65] and of their systematics will addressthese questions and be reported elsewhere. In our studiesof the time evolution of the temperature the consideredsystems did not expand. We plan to investigate the ex-pansion of the system in a next step, see [91], to includeregular contributions into the free energy functional andto extend the treatment of fluctuations to three spatialdimensions. ACKNOWLEDGEMENTS
The authors acknowledge the support of the program“Etoiles montantes en Pays de la Loire 2017”. Thisresearch was supported in part by the ExtreMe Mat-ter Institute (EMMI) at the GSI Helmholtzzentrum f¨urSchwerionenforschung, Darmstadt, Germany. The au-thors thank S. A. Bass and T. Sch¨afer for stimulatingdiscussions.
Appendix A: Numerical implementation
We study the stochastic diffusion equation Eq. (15) bydiscretizing the diffusive net-baryon density on N x sitesequally distributed over a system of longitudinal extent L with resolution (grid spacing) ∆ x = L/N x . Over thefinite ∆ x of a cell, n B must be understood as being av-eraged. Time is discretized in steps of ∆ t at which n B is considered point-wise. Details about the extent anddiscretization in the transverse direction are not impor-tant as we study the evolution of the system and physicalobservables only in the longitudinal direction which is de-coupled from the transverse dynamics. For simplicity weset the transverse area A = 1 fm but have verified theproper behavior with A in the numerics. The Gaussianwhite noise must be understood as averaged over space∆ x and time ∆ t . It is independent between different cellsand time steps with zero mean and variance 1 / (∆ x ∆ t ).Equation (15) is solved by means of a semi-implicitscheme. While the operator associated with the Gaussianmass and surface tension terms is treated implicitly intime, the operator associated with the nonlinear couplingterms is discretized explicitly. The temporal integrationis performed with a predictor-corrector method. For thestochastic diffusion equation of the general form dn B dt = O l n B + O nl ( n B ) + O ξ W (A1)this amounts to solving in a first step (cid:18) − ∆ t O l (cid:19) ˜ n m +1 B = (cid:18) t O l (cid:19) n mB + ∆ t O nl ( n mB )+ ∆ t O ξ W m (A2)for ˜ n m +1 B as an intermediate update of n mB from timestep5 m , and then by using n mB and ˜ n m +1 B in (cid:18) − ∆ t O l (cid:19) n m +1 B = (cid:18) t O l (cid:19) n mB + ∆ t O nl (˜ n mB ) + O nl ( n mB )]+ ∆ t O ξ W m (A3)one finds n m +1 B as the update at timestep m + 1. Theindividual operators in the above equations read (we dropthe index B in the following to improve readability) O l n m = D ∆ x m n c (cid:0) n mj +1 − n mj + n mj − (cid:1) − D ∆ x Kn c (cid:0) n mj +2 − n mj +1 + 6 n mj − n mj − + n mj − (cid:1) , (A4) O nl ( n m ) = D ∆ x (cid:88) i =3 , , λ i n i − c (cid:104) (cid:0) ∆ n mj +1 (cid:1) i − − (cid:0) ∆ n mj (cid:1) i − + (cid:0) ∆ n mj − (cid:1) i − (cid:105) , (A5) O ξ W m = (cid:114) Dn c A ∆ x ∆ t x (cid:16) W mj + − W mj − (cid:17) . (A6)This system of equations is solved for each spatial point j on the grid. The contributions from the nonlinear cou-pling terms are simulated by computing the correspond-ing power of ∆ n mj = n mj − n c at a given site j for eachtimestep m . Without nonlinear coupling terms the pre-dictor and corrector steps are identical and yield exactlythe same solution which makes one of the steps redun-dant. The noise field W in Eq. (A6) has zero mean andvariance 1. Appendix B: Static structure factor in discretizedspace-time
For the Gauss and Gauss+surface models, for whichthe contributions from the nonlinear operator in Eq. (A5)vanish, analytic results for the static structure factor S k in discretized space-time can be derived. In this limit,the general form of the stochastic diffusion equation maybe written in mixed Fourier space as M (1) k ˆ n m +1 k = M ( − k ˆ n mk + N k ˆ W mk (B1)with M ( a ) k = 1 − a ∆ t (cid:16) D ∆ x m n c [2 cos(∆ k ) − − D ∆ x Kn c [2 cos(2∆ k ) − k ) + 6] (cid:17) , (B2) N k = (cid:114) Dn c ∆ tA ∆ x sin(∆ k/ , (B3) ∆ k = k ∆ x and (cid:104) ˆ W mk ( ˆ W mk ) ∗ (cid:105) = 1 /N x . We note thatbased on Eq. (15) we find that Eq. (B1) holds true alsofor the difference ∆ˆ n k instead of ˆ n k . From the definition S mk = V (cid:104) ∆ˆ n mk (∆ˆ n mk ) ∗ (cid:105) (B4)and the condition of stationarity S k = S mk = S m +1 k inequilibrium one finds S k = | N k | A ∆ x (cid:12)(cid:12)(cid:12) M (1) k (cid:12)(cid:12)(cid:12) − (cid:12)(cid:12)(cid:12) M ( − k (cid:12)(cid:12)(cid:12) (B5)= n c m
11 + Km ∆ x [1 − cos(∆ k )] (B6)for the static structure factor, see Eq. (20). This resultis independent of the time-step ∆ t . In the limit of thepure Gaussian model with K = 0 this reduces to S k = n c /m , which is independent of k and ∆ x and agreeswith the result in the continuum, see Eq. (18). We notethat for k = 0 the static structure factor reflects chargeconservation in the entire system. Appendix C: Net-baryon number conservation infinite-size systems
The net-baryon number of the entire system is numer-ically conserved by imposing periodic boundary condi-tions on the net-baryon number density. As a result,charge conservation must be reflected in the behavior ofobservables such as the equal-time correlation function.The latter is connected with the static structure factorby a Fourier transformation. In discretized space onedefines n j = (cid:88) k ˆ n k e ijk ∆ x , (C1)where k = 2 πκ/L is restricted by 0 ≤ κ < N x for κ ∈ Z .Accordingly, the equal-time correlation function followsin discretized space as (cid:104) (∆ n B ) j (∆ n B ) l (cid:105) = 1 V N x − (cid:88) κ =0 e i πκ | j − l | /N x S k . (C2)This definition holds for an infinite system. For the pureGaussian model with S k = n c /m one finds (cid:104) (∆ n B ) j (∆ n B ) l (cid:105) = n c Am δ jl ∆ x (C3)because modes with different κ are orthogonal.For a finite-size system, however, charge conserva-tion must be imposed by demanding that local fluctu-ations vanish upon summation over the entire system,i.e. (cid:80) l (cid:104) (∆ n B ) j (∆ n B ) l (cid:105) = 0 for any j . This is respected6if we impose (cid:104) (∆ n B ) j (∆ n B ) l (cid:105) = 1 V N x − (cid:88) κ =0 e i πκ | j − l | /N x S k − N x V N x − (cid:88) h =0 N x − (cid:88) κ =0 e i πκh/N x S k . (C4)Instead of Eq. (C3) one finds (cid:104) (∆ n B ) j (∆ n B ) l (cid:105) = n c Am (cid:18) δ jl ∆ x − L (cid:19) (C5)for the pure Gaussian model. The finite-size correctionvanishes in the thermodynamic limit for any given reso-lution ∆ x .For the Gauss+surface model with S k given in Eq. (B6)the result of the summations in Eq. (C4) can be obtainednumerically. Due to Eq. (C4) one expects a negativeshift in a finite-size system. This shift has to become lesspronounced with increasing N x , i.e. for fixed resolution∆ x with increasing L . Both these features are seen inthe numerics, cf. Fig. 3. Appendix D: Determination of the correlation length
The form of the equal-time correlation function foundin the equilibrium simulations is that of an exponentialdecay which is modified by a negative shift due to exactcharge conservation. Since in equilibrium the system ishomogeneous on length scales larger than the noise cor-relation this shift is expected to be a constant. Then, asuitable ansatz to determine the correlation length ˜ ξ is (cid:104) (∆ n B ) j (∆ n B ) l (cid:105) = C C exp( −| j − l | ∆ x/C ) + C , (D1)where C is the fit parameter for ˜ ξ in dependence of ∆ x , L and T . The quantity C /C + C gives the value ofthe correlation function over distances of the grid spac-ing ∆ x , i.e. the value of the local variance. The fit resultsfor ˜ ξ shown in Sec. III C are obtained by optimizing thedescription of the local variance in the numerics. We ob-serve that for the Gauss+surface model simulations with L = 20 fm the obtained value of C /C is already quiteclose to the continuum expectation of n c / (2 Am √ K ) forthe local variance in an infinite system, cf. Eqs. (22)and (25), even for T near T c . We note that for smaller L this is not necessarily the case, in particular close to T c .Motivated by the fact that the equilibrium results of theequal-time correlation function in the Ginzburg-Landaumodel can be described by the theoretical expectationof the Gauss+surface model with a modified, effectiveGaussian mass parameter, see Sec. III B, we utilize thesame ansatz Eq. (D1) and strategy in order to fit thenumerical results of the Ginzburg-Landau model and todetermine ˜ ξ . Appendix E: Dynamic structure factor in discretizedspace-time
The diffusion equation in discretized space-timediscussed in Appendix A has for the Gauss andGauss+surface models the following representation in full( ω, k ) Fourier-space: M (1) k e i ∆ ω ˆ n k,ω = M ( − k ˆ n k,ω + N k ˆ W k,ω (E1)with ∆ ω = ω ∆ t . This implies for the correlator (cid:104) ∆ˆ n k,ω ∆ˆ n ∗ k,ω (cid:105) = N k (cid:104) ˆ W k,ω ˆ W ∗ k,ω (cid:105) N ∗ k (cid:16) M (1) k − e − i ∆ ω M ( − k (cid:17) (cid:16) M (1) ∗ k − e i ∆ ω M ( − ∗ k (cid:17) (E2)which gives the dynamic structure factor via S k,ω = lim N t →∞ V ( N t ∆ t ) (cid:104) ∆ˆ n k,ω ∆ˆ n ∗ k,ω (cid:105) , (E3)where N t is the number of (performed) time steps. Fromthis definition it is clear that the dynamic structure factoris a late-time equilibrium observable. For white noise wehave (cid:104) ˆ W k,ω ˆ W ∗ k,ω (cid:105) = 1 / ( N x N t ), and with M ( a ) k and N k defined in Appendix B we obtain S k,ω = 2 n c ˜ χ Dk n c ∆ t − (1 − cos(∆ ω )) + ˜ χ ˜ χ D k (E4)with˜ χ = (1 − cos(∆ k )) / ∆ k , (E5)˜ χ = m (cid:18) − Km ∆ x (cos(∆ k ) − (cid:19) (1 + cos(∆ ω )) . (E6)The result for the pure Gaussian model is found by set-ting K = 0 in Eq. (E6). In the limit of small ∆ ω (cid:28) ω ) and findlim ∆ ω (cid:28) S k,ω = 4 n c ˜ χ Dk ω + 4 D m n c k ˜ χ (cid:0) − Km ∆ x (cos(∆ k ) − (cid:1) . (E7)Moreover, in the limit of small ∆ k (cid:28) χ ≈ − ∆ k and (2 cos(∆ k ) − / ∆ x ≈ − k + k ∆ k in Eq. (E7). Thus, for given ω and k , in the limit of∆ t → x → S ( k, ω ) is recovered from S k,ω . In the numerics, finite resolution in ∆ t and ∆ x implies deviations from the continuum result. Therefore,only the regime of small wavevectors and low frequenciesallows us to judge the accuracy of the numerical scheme.Even in the limit ∆ t →
0, the approach to the continuumis limited to small values of k depending on the spatial7resolution. This limit can be used to determine an ana-lytic expression for the dynamic structure factor S k,t inthe mixed representation from the Fourier transforma-tion into the time-domain. We find S k,t = 12 π (cid:90) ∞−∞ e iωt lim ∆ t → S k,ω dω = S k e − t/τ k , (E8)where S k is the static structure factor in Eq. (B6) and τ k is the relaxation time of fluctuations with wavevector k = 2 πκ/L given via τ − k = 2 Dn c m k ( A k + B k ) (E9) with A k = (1 − cos(∆ k )) / ∆ k , (E10) B k = 2 Kk (1 − cos(∆ k )) m ∆ k . (E11)From Eqs. (E9) - (E11) in the limit of small ∆ x we seethat finite-resolution effects increase τ k compared to thecontinuum result Eq. (28), which is approached in thelimit ∆ x →
0. Moreover, we find that τ k is smaller inthe Gauss+surface model compared to the pure Gaus-sian model with K = 0. This effect is less pronouncedfor small values of κ and away from the transition tem-perature T c . [1] T. Sch¨afer, “Fluid Dynamics and Viscosity in StronglyCorrelated Fluids,” Ann. Rev. Nucl. Part. Sci. , 125-148 (2014) [arXiv:1403.0653 [hep-ph]].[2] P. Kovtun, “Lectures on hydrodynamic fluctuations inrelativistic theories,” J. Phys. A , 473001 (2012)[arXiv:1205.5040 [hep-th]].[3] S. Jeon and U. Heinz, “Introduction to Hydrodynam-ics,” Int. J. Mod. Phys. E , no.10, 1530010 (2015)[arXiv:1503.03931 [hep-ph]].[4] P. Kovtun and L. G. Yaffe, “Hydrodynamic fluctuations,long time tails, and supersymmetry,” Phys. Rev. D ,025007 (2003) [arXiv:hep-th/0303010 [hep-th]].[5] Y. Akamatsu, A. Mazeliauskas and D. Teaney, “A kineticregime of hydrodynamic fluctuations and long time tailsfor a Bjorken expansion,” Phys. Rev. C , no.1, 014909(2017) [arXiv:1606.07742 [nucl-th]].[6] M. Martinez and T. Sch¨afer, “Stochastic hydrodynam-ics and long time tails of an expanding conformalcharged fluid,” Phys. Rev. C , no.5, 054902 (2019)[arXiv:1812.05279 [hep-th]].[7] X. An, G. Basar, M. Stephanov and H. U. Yee, “Rela-tivistic Hydrodynamic Fluctuations,” Phys. Rev. C ,no.2, 024910 (2019) [arXiv:1902.09517 [hep-th]].[8] X. An, G. Basar, M. Stephanov and H. U. Yee, “Fluctua-tion dynamics in a relativistic fluid with a critical point,”[arXiv:1912.13456 [hep-th]].[9] P. C. Hohenberg and B. I. Halperin, “Theory of DynamicCritical Phenomena,” Rev. Mod. Phys. (1977) 435.[10] D. T. Son and M. A. Stephanov, “Dynamic universalityclass of the QCD critical point,” Phys. Rev. D , 056001(2004) [arXiv:hep-ph/0401052 [hep-ph]].[11] H. Fujii and M. Ohtani, “Sigma and hydrodynamicmodes along the critical line,” Phys. Rev. D , 014016(2004) [arXiv:hep-ph/0402263 [hep-ph]].[12] B. V. Jacak and B. M¨uller, “The exploration of hot nu-clear matter,” Science , 310-314 (2012).[13] P. Braun-Munzinger, V. Koch, T. Sch¨afer and J. Stachel,“Properties of hot and dense matter from relativisticheavy ion collisions,” Phys. Rept. , 76-126 (2016)[arXiv:1510.00442 [nucl-th]].[14] W. Busza, K. Rajagopal and W. van der Schee, “HeavyIon Collisions: The Big Picture, and the Big Ques-tions,” Ann. Rev. Nucl. Part. Sci. , 339-376 (2018)[arXiv:1802.04801 [hep-ph]]. [15] D. A. Teaney, “Viscous Hydrodynamics and the QuarkGluon Plasma,” [arXiv:0905.2433 [nucl-th]].[16] B. Schenke, S. Jeon and C. Gale, “(3+1)D hydrodynamicsimulation of relativistic heavy-ion collisions,” Phys. Rev.C , 014903 (2010) [arXiv:1004.1408 [hep-ph]].[17] U. Heinz and R. Snellings, “Collective flow and viscosityin relativistic heavy-ion collisions,” Ann. Rev. Nucl. Part.Sci. , 123-151 (2013) [arXiv:1301.2826 [nucl-th]].[18] L. Del Zanna, V. Chandra, G. Inghirami, V. Rolando,A. Beraudo, A. De Pace, G. Pagliara, A. Drago andF. Becattini, “Relativistic viscous hydrodynamics forheavy-ion collisions with ECHO-QGP,” Eur. Phys. J. C , 2524 (2013) [arXiv:1305.7052 [nucl-th]].[19] I. Karpenko, P. Huovinen and M. Bleicher, “A 3+1dimensional viscous hydrodynamic code for relativisticheavy ion collisions,” Comput. Phys. Commun. ,3016-3027 (2014) [arXiv:1312.4160 [nucl-th]].[20] R. Derradi de Souza, T. Koide and T. Kodama, “Hy-drodynamic Approaches in Relativistic Heavy Ion Re-actions,” Prog. Part. Nucl. Phys. , 35-85 (2016)[arXiv:1506.03863 [nucl-th]].[21] P. Romatschke and U. Romatschke, “Relativistic FluidDynamics In and Out of Equilibrium,” [arXiv:1712.05815[nucl-th]].[22] M. Connors, C. Nattrass, R. Reed and S. Salur, “Jetmeasurements in heavy ion physics,” Rev. Mod. Phys. , 025005 (2018) [arXiv:1705.01974 [nucl-ex]].[23] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Sz-abo, “The Order of the quantum chromodynamics transi-tion predicted by the standard model of particle physics,”Nature , 675-678 (2006) [arXiv:hep-lat/0611014 [hep-lat]].[24] M. M. Aggarwal et al. [STAR], “An Experimental Ex-ploration of the QCD Phase Diagram: The Search forthe Critical Point and the Onset of De-confinement,”[arXiv:1007.2613 [nucl-ex]].[25] B. Friman, C. H¨ohne, J. Knoll, S. Leupold, J. Ran-drup, R. Rapp and P. Senger, “The CBM physics book:Compressed baryonic matter in laboratory experiments,”Lect. Notes Phys. , pp.1-980 (2011).[26] X. Luo, “Exploring the QCD Phase Structure with BeamEnergy Scan in Heavy-ion Collisions,” Nucl. Phys. A ,75-82 (2016) [arXiv:1512.09215 [nucl-ex]]. [27] A. Bzdak, S. Esumi, V. Koch, J. Liao, M. Stephanovand N. Xu, “Mapping the Phases of Quantum Chromo-dynamics with Beam Energy Scan,” Phys. Rept. ,1-87 (2020) [arXiv:1906.00936 [nucl-th]].[28] X. Luo, S. Shi, N. Xu and Y. Zhang, “A Study of theProperties of the QCD Phase Diagram in High-EnergyNuclear Collisions,” Particles , no.2, 278-307 (2020)[arXiv:2004.00789 [nucl-ex]].[29] K. Rajagopal and F. Wilczek, “Static and dynamic crit-ical phenomena at a second order QCD phase transi-tion,” Nucl. Phys. B , 395-425 (1993) [arXiv:hep-ph/9210253 [hep-ph]].[30] J. Berges and K. Rajagopal, “Color superconductiv-ity and chiral symmetry restoration at nonzero baryondensity and temperature,” Nucl. Phys. B , 215-232(1999) [arXiv:hep-ph/9804233 [hep-ph]].[31] A. M. Halasz, A. D. Jackson, R. E. Shrock,M. A. Stephanov and J. J. M. Verbaarschot, “On thephase diagram of QCD,” Phys. Rev. D , 096007 (1998)[arXiv:hep-ph/9804290 [hep-ph]].[32] K. Fukushima and T. Hatsuda, “The phase diagramof dense QCD,” Rept. Prog. Phys. (2011), 014001[arXiv:1005.4814 [hep-ph]].[33] K. Fukushima and C. Sasaki, “The phase diagram of nu-clear and quark matter at high baryon density,” Prog.Part. Nucl. Phys. (2013), 99-154 [arXiv:1301.6377[hep-ph]].[34] M. A. Stephanov, K. Rajagopal and E. V. Shuryak, “Sig-natures of the tricritical point in QCD,” Phys. Rev. Lett. , 4816-4819 (1998) [arXiv:hep-ph/9806219 [hep-ph]].[35] M. A. Stephanov, K. Rajagopal and E. V. Shuryak,“Event-by-event fluctuations in heavy ion collisions andthe QCD critical point,” Phys. Rev. D , 114028 (1999)[arXiv:hep-ph/9903292 [hep-ph]].[36] Y. Hatta and M. A. Stephanov, “Proton number fluctua-tion as a signal of the QCD critical endpoint,” Phys. Rev.Lett. , 102003 (2003) [arXiv:hep-ph/0302002 [hep-ph]].[37] M. Asakawa and M. Kitazawa, “Fluctuations of con-served charges in relativistic heavy ion collisions: An in-troduction,” Prog. Part. Nucl. Phys. , 299-342 (2016)[arXiv:1512.05038 [nucl-th]].[38] X. Luo and N. Xu, “Search for the QCD Critical Pointwith Fluctuations of Conserved Quantities in RelativisticHeavy-Ion Collisions at RHIC : An Overview,” Nucl. Sci.Tech. , no.8, 112 (2017) [arXiv:1701.02105 [nucl-ex]].[39] L. Adamczyk et al. [STAR], “Energy Dependenceof Moments of Net-proton Multiplicity Distributionsat RHIC,” Phys. Rev. Lett. , 032302 (2014)[arXiv:1309.5681 [nucl-ex]].[40] J. Adam et al. [STAR], “Net-proton number fluctua-tions and the Quantum Chromodynamics critical point,”[arXiv:2001.02852 [nucl-ex]].[41] J. Adamczewski-Musch et al. [HADES], “Proton num-ber fluctuations in √ s NN = 2.4 GeV Au+Au collisionsstudied with HADES,” [arXiv:2002.08701 [nucl-ex]].[42] M. A. Stephanov, “Non-Gaussian fluctuations near theQCD critical point,” Phys. Rev. Lett. , 032301 (2009)[arXiv:0809.3450 [hep-ph]].[43] M. Asakawa, S. Ejiri and M. Kitazawa, “Third momentsof conserved charges as probes of QCD phase structure,”Phys. Rev. Lett. , 262301 (2009) [arXiv:0904.2089[nucl-th]].[44] M. A. Stephanov, “On the sign of kurtosis near theQCD critical point,” Phys. Rev. Lett. , 052301 (2011) [arXiv:1104.1627 [hep-ph]].[45] B. Berdnikov and K. Rajagopal, “Slowing out-of-equilibrium near the QCD critical point,” Phys. Rev. D , 105017 (2000) [arXiv:hep-ph/9912274 [hep-ph]].[46] M. Nahrgang, S. Leupold, C. Herold and M. Bleicher,“Nonequilibrium chiral fluid dynamics including dissi-pation and noise,” Phys. Rev. C , 024912 (2011)[arXiv:1105.0622 [nucl-th]].[47] M. Nahrgang, S. Leupold and M. Bleicher, “Equilibra-tion and relaxation times at the chiral phase transitionincluding reheating,” Phys. Lett. B , 109-116 (2012)[arXiv:1105.1396 [nucl-th]].[48] M. Nahrgang, C. Herold, S. Leupold, I. Mishustin andM. Bleicher, “The impact of dissipation and noise on fluc-tuations in chiral fluid dynamics,” J. Phys. G , 055108(2013) [arXiv:1105.1962 [nucl-th]].[49] C. Herold, M. Nahrgang, I. Mishustin and M. Bleicher,“Chiral fluid dynamics with explicit propagation of thePolyakov loop,” Phys. Rev. C , no.1, 014907 (2013)[arXiv:1301.1214 [nucl-th]].[50] M. Nahrgang, C. Herold and M. Bleicher, “Influence of aninhomogeneous and expanding medium on signals of theQCD phase transition,” Nucl. Phys. A , 899c-902c (2013) [arXiv:1301.2577 [nucl-th]].[51] C. Herold, M. Nahrgang, Y. Yan and C. Kobdaj, “Net-baryon number variance and kurtosis within nonequilib-rium chiral fluid dynamics,” J. Phys. G , no.11, 115106(2014) [arXiv:1407.8277 [hep-ph]].[52] S. Mukherjee, R. Venugopalan and Y. Yin, “Real timeevolution of non-Gaussian cumulants in the QCD crit-ical regime,” Phys. Rev. C , no.3, 034912 (2015)[arXiv:1506.00645 [hep-ph]].[53] C. Herold, M. Nahrgang, Y. Yan and C. Kobdaj,“Dynamical net-proton fluctuations near a QCD crit-ical point,” Phys. Rev. C , no.2, 021902 (2016)[arXiv:1601.04839 [hep-ph]].[54] M. Nahrgang and C. Herold, “Phenomena at the QCDphase transition in nonequilibrium chiral fluid dynam-ics (N χ FD),” Eur. Phys. J. A , no.8, 240 (2016)[arXiv:1602.07223 [nucl-th]].[55] S. Mukherjee, R. Venugopalan and Y. Yin, “Universaloff-equilibrium scaling of critical cumulants in the QCDphase diagram,” Phys. Rev. Lett. , no.22, 222301(2016) [arXiv:1605.09341 [hep-ph]].[56] C. Herold, M. Bleicher, M. Nahrgang, J. Steinheimer,A. Limphirat, C. Kobdaj and Y. Yan, “Broadeningof the chiral critical region in a hydrodynamically ex-panding medium,” Eur. Phys. J. A , no.2, 19 (2018)[arXiv:1710.03118 [hep-ph]].[57] M. Stephanov and Y. Yin, “Hydrodynamics with para-metric slowing down and fluctuations near the crit-ical point,” Phys. Rev. D , no.3, 036006 (2018)[arXiv:1712.10305 [nucl-th]].[58] C. Herold, A. Kittiratpattana, C. Kobdaj, A. Limphi-rat, Y. Yan, M. Nahrgang, J. Steinheimer and M. Ble-icher, “Entropy production and reheating at the chiralphase transition,” Phys. Lett. B , 557-562 (2019)[arXiv:1810.02504 [hep-ph]].[59] K. Rajagopal, G. Ridgway, R. Weller and Y. Yin, “Hy-dro+ in Action: Understanding the Out-of-EquilibriumDynamics Near a Critical Point in the QCD Phase Dia-gram,” [arXiv:1908.08539 [hep-ph]].[60] L. Du, U. Heinz, K. Rajagopal and Y. Yin, “Fluc-tuation dynamics near the QCD critical point,” [arXiv:2004.02719 [nucl-th]].[61] M. Kitazawa, M. Asakawa and H. Ono, “Non-equilibriumtime evolution of higher order cumulants of conservedcharges and event-by-event analysis,” Phys. Lett. B ,386-392 (2014) [arXiv:1307.2978 [nucl-th]].[62] M. Sakaida, M. Asakawa and M. Kitazawa, “Effectsof global charge conservation on time evolution of cu-mulants of conserved charges in relativistic heavy ioncollisions,” Phys. Rev. C , no.6, 064911 (2014)[arXiv:1409.6866 [nucl-th]].[63] M. Sakaida, M. Asakawa, H. Fujii and M. Kitazawa, “Dy-namical evolution of critical fluctuations and its obser-vation in heavy ion collisions,” Phys. Rev. C , no.6,064905 (2017) [arXiv:1703.08008 [nucl-th]].[64] M. Nahrgang, M. Bluhm, T. Sch¨afer and S. A. Bass,“Baryon number diffusion with critical fluctuations,”Nucl. Phys. A , 824-827 (2017) [arXiv:1804.02976[nucl-th]].[65] M. Nahrgang, M. Bluhm, T. Sch¨afer and S. A. Bass,“Diffusive dynamics of critical fluctuations near the QCDcritical point,” Phys. Rev. D , no.11, 116015 (2019)[arXiv:1804.05728 [nucl-th]].[66] M. Bluhm, Y. Jiang, M. Nahrgang, J. Pawlowski, F. Ren-necke and N. Wink, “Time-evolution of fluctuations assignal of the phase transition dynamics in a QCD-assistedtransport approach,” Nucl. Phys. A , 871-874 (2019)[arXiv:1808.01377 [hep-ph]].[67] Y. Akamatsu, D. Teaney, F. Yan and Y. Yin, “Transits ofthe QCD critical point,” Phys. Rev. C , no.4, 044901(2019) [arXiv:1811.05081 [nucl-th]].[68] M. Bluhm and M. Nahrgang, “Time-evolution of net-baryon density fluctuations across the QCD critical re-gion,” [arXiv:1911.08911 [nucl-th]].[69] M. Bluhm, M. Nahrgang, A. Kalweit, M. Arslan-dok, P. Braun-Munzinger, S. Floerchinger, E. S. Fraga,M. Gazdzicki, C. Hartnack, C. Herold, R. Holzmann,I. Karpenko, M. Kitazawa, V. Koch, S. Leupold,A. Mazeliauskas, B. Mohanty, A. Ohlson, D. Oliiny-chenko, J. M. Pawlowski, C. Plumberg, G. W. Ridgway,T. Sch¨afer, I. Selyuzhenkov, J. Stachel, M. Stephanov,D. Teaney, N. Touroux, V. Vovchenko and N. Wink, “Dy-namics of critical fluctuations: Theory phenomenologyheavy-ion collisions,” [arXiv:2001.08831 [nucl-th]].[70] C. Nonaka and M. Asakawa, “Hydrodynamical evolu-tion near the QCD critical end point,” Phys. Rev. C ,044904 (2005) [arXiv:nucl-th/0410078 [nucl-th]].[71] M. Bluhm and B. K¨ampfer, “Quasi-particle perspec-tive on critical end-point,” PoS CPOD2006 , 004 (2006)[arXiv:hep-ph/0611083 [hep-ph]].[72] P. Parotto, M. Bluhm, D. Mroczek, M. Nahrgang,J. Noronha-Hostler, K. Rajagopal, C. Ratti, T. Sch¨aferand M. Stephanov, “QCD equation of state matched tolattice data and exhibiting a critical point singularity,”Phys. Rev. C , no.3, 034901 (2020) [arXiv:1805.05249[hep-ph]].[73] C. Young, “Numerical integration of thermal noise in rel-ativistic hydrodynamics,” Phys. Rev. C (2014) no.2,024913 [arXiv:1306.0472 [nucl-th]].[74] K. Murase and T. Hirano, “Hydrodynamic fluctuationsand dissipation in an integrated dynamical model,” Nucl.Phys. A , 276-279 (2016) [arXiv:1601.02260 [nucl-th]].[75] M. Nahrgang, M. Bluhm, T. Sch¨afer and S. Bass, “To-ward the description of fluid dynamical fluctuations inheavy-ion collisions,” Acta Phys. Polon. Supp. , 687 (2017) [arXiv:1704.03553 [nucl-th]].[76] M. Bluhm, M. Nahrgang, T. Sch¨afer and S. A. Bass,“Fluctuating fluid dynamics for the QGP in the LHCand BES era,” EPJ Web Conf. , 16004 (2018)[arXiv:1804.03493 [nucl-th]].[77] M. Singh, C. Shen, S. McDonald, S. Jeon and C. Gale,“Hydrodynamic Fluctuations in Relativistic Heavy-Ion Collisions,” Nucl. Phys. A , 319-322 (2019)[arXiv:1807.05451 [nucl-th]].[78] T. Hirano, R. Kurita and K. Murase, “Hydrody-namic fluctuations of entropy in one-dimensionally ex-panding system,” Nucl. Phys. A , 44-67 (2019)[arXiv:1809.04773 [nucl-th]].[79] A. Sakai, K. Murase and T. Hirano, “Rapidity decorre-lation of anisotropic flow caused by hydrodynamic fluc-tuations,” [arXiv:2003.13496 [nucl-th]].[80] J. B. Bell, A. L. Garcia and S. A. Williams, “Numeri-cal Methods for the Stochastic Landau-Lifshitz Navier-Stokes Equations,” Phys. Rev. E , 016708 (2007)[arXiv:math/0612324 [math.NA]].[81] A. Donev, E. Vanden-Eijnden, A. L. Garcia andJ. B. Bell, “On the Accuracy of Finite-Volume Schemesfor Fluctuating Hydrodynamics,” [arXiv:0906.2425[physics.flu-dyn]].[82] J. A. de la Torre, P. Espa˜nol and A. Donev, “Finite ele-ment discretization of non-linear diffusion equations withthermal fluctuations,” J. Chem. Phys. (2015) 094115[arXiv:1410.6340 [cond-mat.stat-mech]].[83] M. M. Tsypin, “Universal effective potential for scalarfield theory in three-dimensions by Monte Carlo compu-tation,” Phys. Rev. Lett. , 2015-2018 (1994).[84] M. M. Tsypin, “Effective potential for a scalar fieldin three dimensions: Ising model in the ferromagneticphase,” Phys. Rev. B , 8911-8917 (1997).[85] R. Guida and J. Zinn-Justin, “3-D Ising model: The Scal-ing equation of state,” Nucl. Phys. B , 626-652 (1997)[arXiv:hep-th/9610223 [hep-th]].[86] M. Agah Nouhou, M. Bluhm, A. Borer, M. Nahrgang,T. Sami and N. Touroux, “Finite size effects on cu-mulants of the critical mode,” PoS CORFU2018 , 179(2019) [arXiv:1906.02647 [nucl-th]].[87] M. Bluhm, M. Nahrgang, S. A. Bass and T. Sch¨afer,“Impact of resonance decays on critical point signals innet-proton fluctuations,” Eur. Phys. J. C , no.4, 210(2017) [arXiv:1612.03889 [nucl-th]].[88] M. Bluhm, M. Nahrgang, S. A. Bass and T. Sch¨afer,“Behavior of universal critical parameters in the QCDphase diagram,” J. Phys. Conf. Ser. , no.1, 012074(2017) [arXiv:1612.04564 [nucl-th]].[89] J. Randrup and J. Cleymans, “Exploring high-densitybaryonic matter: Maximum freeze-out density,” Eur.Phys. J. , 218-219 (2016) [arXiv:0905.2824 [nucl-th]].[90] L. V. Bravina, I. Arsene, M. S. Nilsson, K. Tywo-niuk, E. E. Zabrodin, J. Bleibel, A. Faessler, C. Fuchs,M. Bleicher, G. Burau and H. St¨ocker, “Microscopicmodels and effective equation of state in nuclear colli-sions at FAIR energies,” Phys. Rev. C78