The equilibrium and dynamical cumulants of QCD chiral order parameter with parametric Landau free energy

By linearly parameterizing the QCD Landau free energy near the critical point in the baryon chemical potential and temperature plane, we study the fluctuations of the QCD chiral order parameter field (the σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} field) in the equilibrium case and dynamical phase transition, respectively. By setting the system size to the typical size of the QGP fireball (≈103\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx 10^3$$\end{document} fm3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^3$$\end{document}), we show that in the equilibrium case, the discontinuity of the order parameter in the first order phase transition region is replaced by smooth crossover, and the corresponding fluctuations are broadened. Meanwhile, the quartic cumulant κ4\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa _4$$\end{document} of the σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} field is generally negative near the phase transition line. We further derive the dynamical evolution of the QCD Landau free energy in the Fokker–Plank framework, based on which we deduce the dynamical cumulants of the σ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sigma $$\end{document} field. Assuming the temperature decreases as a known function of time, we numerically evaluate the dynamical cumulants and confirm that the cumulants present clear memory effects. Moreover, the memory effects on the first order phase transition side is stronger than that on the crossover side, and the dynamical cumulants at the hypothetical freeze-out line present rich non-monotonic structures.


Introduction
The QCD phase structure has been of intensive theoretical and experimental interest for decades [1,2,3,4,5,6]. In various extreme conditions, theoretical studies predict different phases of the QCD matter, such as quark-gluon plasma (QGP), hadronic resonance gas (HRG), and color superconductor. Of great interest is the chiral phase transition (χPT) between the hot QGP phase and the HRG phase [7,9,8]. For small baryon chemical potential µ, lattice simulations show that the transition is a broad crossover [10,11,12,13]. For large µ, the sign problem of lattice QCD prevents full ab initio simulations. The phenomenological models such as the Nambu-Jona-Lasinio model [14,15,16,17] and quark-meson model [18,19,20], and nonperturbative methods like Functional Renormalization Group method [21] and Dyson-Schwinger equation [22,23] provide relatively complete description on the χPT, predicting crossover at small µ, first order phase transition at large µ, and the existence of a QCD critical point.
On the other hand, owing to the strong coupling between the chiral σ field and quarks in QCD-inspired models, the highorder cumulants of net-proton production are expected to be sensitive to the increase of fluctuations of the σ field near the critical point [24,25]. The related observables are being measured by the STAR collaboration at the Relativistic Heavy Ion Collider (RHIC) [26,27,28]. The latest experimental data of κσ 2 for the net-proton production presents a non-monotonic a e-mail: junhui.zheng@nwu.edu.cn (corresponding author). variation as a function of collision energy for the region √ s NN = 7.7 − 200 GeV in the Au+Au central collisions [28], which partially agrees with theoretical expectations [24]. However, as shown in Ref. [29], by introducing a freeze-out scheme to the hydrodynamics, the model calculations of the equilibrium critical fluctuations of the net protons still failed to qualitatively explain the full experimental data. It presented the limitation of thermal equilibrium assumption and suggested that the critical dynamics should be taken into account for the description of dynamical χPT in heavy ion collision.
In recent years, to study the dynamics of non-equilibrium fluctuations, including the dynamical evolution of the order parameter field and the diffusion of the conserved charges, different dynamical models have been developed [30,31,32,33,34,35,36]. The dynamical critical fluctuations consistently present clear memory effects and critical slowing down effects. As a result, both the sign and the magnitude of the high order cumulants can be different from the equilibrium ones, and the kurtosis shows non-monotonic behaviors. To be closer to the experimental process in RHIC, the coupled dynamical evolution of the critical modes and the hydrodynamic background are further developed, like that has been done in the chiral hydrodynamics [37,38,39,40,41] and hydro+ [42,43,44,45]. Moreover, other kinds of factors such as the spatially nonuniform temperature (and chemical potential) effects on the fluctuations of the order parameter [46], the non-critical fluctuations [47], and the proper freeze-out scheme [29,48] will also influence the theoretical predictions significantly.
As a main component of the critical dynamics, parameterization of the QCD equation of state through the Ising mapping is usually applied to get insight into the universal behaviors in the critical region [49,50,51,52]. However, due to the discontinuities of QCD equation of state on the first order phase transition side, the critical dynamics remains unclear and requires more exploration. In this paper, we focus on the dynamics of critical fluctuations in the χPT region (including the first order phase transition and crossover) rather than building the complete dynamical modeling. We provide an alternative way to parameterize the QCD Landau free energy in different phase transition scenarios (which also serves as the basis of our discussion on the spatially nonuniform temperature effects in Ising-like models [46]), evaluate the finite size effects of the QCD matter on the order parameter and fluctuations, and develop a distinct set of dynamical equations based on the Fokker-Plank equation to study the dynamical free energy and the dynamical cumulants. The article is structured as following. In Sec. 2, we linearly parameterize the Landau free energy in the (µ, T ) plane, where both the crossover and the first order phase transition side are described. In Sec. 3, we set up the parameters in the free energy and present the results of the equilibrium cumulants in the system of different volumes. In Sec. 4, we deduce the dynamical free energy based on the Fokker-Plank equation. In Sec. 5, we numerically calculate the non-equilibrium cumulants at different scenarios and on a hypothetical freeze-out line with a fixed volume V = 10 3 fm 3 . Finally, in Sec. 6, we summarize the main results of this paper and give a discussion.

Parameterization of the free energy
Because the critical behaviors of the σ field fall into the same universality class as the 3D Ising model [53], in principle, one can parameterize the QCD equation of state by directly mapping the QCD parameters (µ, T ) to the Ising variables (r, h), as is done in Refs. [54,55,56,57]. Combined with the dynamical models, there are various interesting dynamical effects presented, like memory effects and universal off-equilibrium scaling behaviors [55,31,58,36]. However, the early studies focus their discussions of critical dynamics on the crossover side, due to the discontinuities of the equations of state in the first order phase transition region. In this section, we develop an alternative method to directly parameterize the QCD free energy, in which the two phase transition scenarios are unified in a same framework.
According to the Landau theory of phase transition, the free energy in the critical region is supposed to be analytic and obeys the symmetry of the Hamiltonian. Then, the Landau free energy density of the χPT [59,60] can be generally written in terms of the σ field as which is Taylor expanded as a function of σ up to the fourth order. The constant term is omitted because it is irrelevant to the structure of the QCD phase diagram. Note that a zero-momentum mode approximation of the σ field is assumed in above, be- cause we focus on the phase transition region in the vicinity of the critical point, where the long-wavelength modes are dominant. The distribution of the critical field is described by the probability distribution function P[σ] ∝ exp{−Ω[σ]V/T } [24], where V is the volume of the system. In the chiral limit, α 1 = α 3 = 0, while for the physical world, a finite α 1 σ term is introduced to handle the explicit chiral symmetry breaking of the quark masses [61]. The cubic term, α 3 σ 3 , emerges after the renormalization contributed from the high-momentum modes of the σ field. The coefficient of the quartic term, α 4 , is supposed to be positive, in order to sustain the stability of the system.
In the thermodynamic limit (V → ∞), the phase structure is fully determined by the global minimum of the free energy (1). It is convenient to present the information of the χPT by performing a translation transformation for the σ field in the free energy, i.e, σ =σ + σ c , where Consequently, the cubic term in Eq. (1) can be eliminated, and we obtain the translated free energy density (after getting rid of a constant term) which has exactly the same form as the free energy density of the Ising model. The redefined coefficients, η i , can be expressed as functions of the original coefficients α i and the shift σ c by Now theσ field has become an effective order parameter, which behaves just like the spin density in the Ising model. With the translated free energy, we can divide the phase spaces in the µ-T plane into four zones, according barely to the sign of η 1 and η 2 . As shown in Fig. 1, the phase transition line is determined by η 1 = 0, and the order of the phase transition is determined by the sign of η 2 . At the right part of the diagram (zone II and IV), η 2 < 0, the phase transition is first order, as there is a coexistence region around the phase transition line, in which two local minima of the free energy coexist. The coexistence region is separated out by the dash-dotted lines, which are determined by ∆ ≡ ( η 1 2η 4 ) 2 + ( η 2 3η 4 ) 3 < 0. The physical vacuum locates at the lower minimum of the free energy. This is the typical feature of the first order phase transition. For the other parts of the µ-T plane where ∆ > 0, only one minimum exists, and the system undergoes a continuous phase transition (crossover) as η 1 crosses its zero point. The sign of η 1 determines the location of the global minimum. In general, in the upper regions (zone I and II), with η 1 > 0, the global minimum is atσ < 0 (i.e. σ < σ c ); while in the lower regions (zone III and IV), with η 1 < 0, the global minimum is atσ > 0 (i.e. σ > σ c ). At the critical point (µ c , T c ), both η 1 and η 2 vanish.
To parameterize the coefficients in Ω[σ] linearly, it is convenient to define two unit vectors b = (cos θ b , sin θ b ) and b = (cos θ b , sin θ b ) in the µ-T plane. The vector b is parallel to the tangent of the phase transition curve at the critical point; The vector b is along the boundary of region I and II in Fig. 1. Note that b and b are not necessarily orthogonal. The angle between them are determined by the realistic QCD equation of state around the critical point, which is still under study. By projecting the vector (µ−µ c , T −T c ) onto the perpendicular vector of b and that of b , respectively, the coefficients are linearly parameterized as: With these projections, the sign of η 1 and η 2 at different parts (I-IV) can be correctly expressed by simply constraining all the constants to d i > 0 (i = 1, 2, 4), i.e., in the right side of the vector b, η 1 > 0, and in the left side of b , η 2 > 0. Note that the magnitude of d i are again determined by the QCD equation of state, which currently remains unknown. In the following, for simplicity and illustration, we treat them as input parameters, and make qualitative calculations and discussions. Along the phase transition curve, η 1 (µ, T ) = 0, the correlation length is ξ = η −1/2 2 on the crossover side and ξ = (−2η 2 ) −1/2 on the first order phase transition side, which diverges at the critical point where η 2 = 0 [4]. In the critical region, we assume that the change of σ c on the phase transition line is small and treat the variable σ c (µ, T ) as a constant, σ c (µ, T ) = σ c (µ c , T c ), in the zero-order approximation of variables µ − µ c and T − T c . Then, the translated free energy density (3) is fully parameterized by µ c , T c , σ c , d 1 , d 2 , d 4 , θ b , and θ b . By using the relations Eq. (2) and Eqs. (4)-(6), the original free energy Ω[σ] (Eq.(1)) can also be expressed in the parameterized form.
Through the transformation from Eq. (1) to Eq. (3) and the parameterization Eqs. (7)-(9), we have built the connection between the QCD free energy and the Ising free energy at the mean-field level. Indeed, Eqs. (7) and (8) describe the dependence of η 1 and η 2 on the QCD variables (µ, T ) in the linear approximation, and the coefficients η 1 and η 2 in the translated free energy (Eq. (3)) are directly related to the reduced magnetic field and temperature h and r in the Ising model, respectively. Note that the current parameterization of the QCD free energy is only up to the linear order, which works in the region close to the critical point. When treating a larger region in the µ-T plane, higher order expansions of µ − µ c , and T − T c should be taken into account. It is worth stressing that this parameterization can also be applied to the (hadronic) gas-liquid phase transition [62,63] or to other similar scenarios, as long as they belong to the same universality class. With this parametric free energy, we can easily develop approaches to study the equilibrium and non-equilibrium cumulants in different phase transition scenarios, as will be shown in the following sections.

parameter setup and equilibrium cumulants
In this section, we calculate the equilibrium cumulants based on the above parametric Landau free energy. The cumulants of the σ field are defined as where the moments µ n are and the probability distribution function for the σ field is The numerical calculation of the cumualnts are implemented with the following setup of parameters: (µ c , T c ) = (240, 170) MeV, σ c = 50 MeV, d 1 = 3 × 10 4 MeV 2 , d 2 = 400 MeV, d 4 = 15, sin θ b = − cos θ b = 0.99, and cos θ b = sin θ b = 0.141. Here, b ⊥ b is set for simplicity. The choice of these parameters are constrained by comparing with the effective potential and cumulants from the linear sigma model with constituent quarks in [37], with which the equilibrium values of the cumulants in the phase transition region are approximately of the order κ n ∼ 10 n MeV n .
Note that the discussion and parameterizations of the free energy density (1) in above is in the thermodynamic limit. In reality, especially in the experimental environment, the size of the QGP fireball in the heavy ion collision is finite (the typical volume of the fireball is approximately V = 10 3 fm 3 ). The finite size has a direct influence on the renormalization process. Consequently, the coefficients α i in Eq. (1) depend on the volume, and further, the parameters µ c , T c , σ c , d 1 , d 2 , d 4 , θ b , and θ b are also volume-dependent. There is a standard process for the renormalization of the coefficients of the free energy [64]. However, even through the volume will change the magnitude of the parameters, the parameterized form of the Landau free energy density will not be changed. Therefore, in this article, we will not discuss the finite size corrections to the free energy density itself but simply utilize the same parameter set for different volumes to study the behaviors of the cumulants. Even without considering the influence on the free energy density from the renormalization of size, we still have to face the finite-volume effects on the fluctuations of the σ field as shown in the probability function (15). In Fig. 2(a), we present the first two order cumulants of the σ field with respect to temperature for different volumes. The chemical potential is fixed at µ = 250 MeV so the system undergoes a first order phase transition as the change of temperature. For a large enough volume (say 10 5 fm 3 ), the expectation value κ 1 = σ is discontinuous because the barrier of the free energy (∆ΩV) can not be overcome by the thermodynamic fluctuations, i.e., ∆ΩV T . The σ field locates only near the global minimum point and the variance κ 2 is also suppressed by the volume. Specifically, on the phase transition line η 1 = 0, the correlation length of the σ field is ξ = (−2η 2 ) −1/2 and the barrier is ∆ΩV = η 2 2 V/4η 4 . Thus, the discontinuity of the first order phase transition becomes evident when V 16T η 4 ξ 4 . For a smaller system volume, the value of T/V in the probability distribution function P[σ] grows. As a result, in the first order phase transition region, the σ field is allowed to be at both of the local minimum points of the free energy with considerable probability, which leads to the enhancements of variance κ 2 . Meanwhile, the expectation value κ 1 deviates from the global minimum of the free energy, and its discontinuity is rounded. Similar finite-size rounding effects of a first order phase transition can also be found in Refs. [65,66,67]. This hints that for a small system, the perturbation theory around the global minimum does not work well anymore, and the equation of states obtained in the thermodynamic limit [24,31] are not suitable for the system with small size.
In Fig. 2(b), we also plot κ 2 -κ 4 as functions of the volume of the system for (T, µ) locating in the crossover regime, and make comparisons with the Gaussian fluctuationsκ 2 = T ξ 2 /V and the perturbative non-Gaussian cumulantsκ 3 = −2λ 3 (T/V) 2 ξ 6 andκ 4 = 6 (T/V) 3 [2(λ 3 ξ) 2 − λ 4 ]ξ 8 [31]. For large volume, the variance κ 2 is consistent withκ 2 . But as the decreasing of the volume, the non-Gaussian corrections are enhanced as the widening of the variance, and lead to a significant deviation of κ 2 fromκ 2 . The Gaussian approximation is valid when the non-Gaussian contributions are small. Expanding Eq. (3) around the extreme pointσ m where ∂σΩ = 0, the third and forth order coupling constants of the fluctuation δσ defined as in Refs. [24,31] become λ 3 = 3η 4σm and λ 4 = η 4 . The corresponding conditions for the validity of the Gaussian approximation are λ 3 (δσ) 3 V/3T 1 and λ 4 (δσ) 4 V/4T 1. The typical magnitude of δσ is of the order of √κ 2 . Thus, the conditions become V (η 4σm ) 2 ξ 6 T and V η 4 ξ 4 T/4. Specially on the phase transition line η 1 = 0, we haveσ m = 0 and the first condition is satisfied automatically. In the condition of large volume, the higher order cumulants κ 3 and κ 4 also become consistent with the perturbative results as shown in the figure.
In the following calculations, we fix the system volume to the typical size of the QGP fireball, V = 10 3 fm 3 . In Fig. 3, we present the density plots of the equilibrium cumulants κ 2 − κ 4 in the η 1 -η 2 plane. The magnitude of κ 2 , κ 3 , and κ 4 in the first order phase transition region (η 2 < 0) are generally larger than that in the crossover region (η 2 > 0) for a given small η 1 . Similar to Ref. [24], κ 3 changes its sign after crossing the phase transition line. On the other hand, unlike the results in Ref. [24], κ 4 is generally negative (red color region) near the phase transition line (η 1 = 0), for both the crossover side and the first order phase transition side. The change of sign for κ 4 happens at the dotted blue lines, inside which the red color region is much larger than the region rounded up by the dashed purple lines (refers to the sign-change line of κ 4 in the thermodynamic limit). As explained in the above paragraphs, the broadening of the negative κ 4 region on the first order phase transition side is due to the rounding effects, since the twopeak shape of P[σ] (for the double-well of Ω[σ] with small volume) has less kurtosis than the Gaussian distribution [57]. While in the thermodynamic limit for the first order phase transition region, even through the free energy Ω[σ] also has the double-well shape, one peak of P[σ] is strongly enhanced by the large volume comparing to the other one, thus P[σ] shapes like a one-peak distribution and only the region near the global minimum of Ω[σ] is distributed. As a result, κ 4 is negative only in the crossover region in the thermodynamic limit [57], striking contrast to the results in systems of small size.

Dynamical equations for the free energy
In this section, we employ the Fokker-Plank equation to study the dynamical behaviors of the parameterized QCD free en- Fig. 3: The equilibrium cumulants κ 2 − κ 4 in the η 1 -η 2 plane for V = 10 3 fm 3 . The colors red and blue refer to negative and positive sign of the cumulants, respectively, and darker colors correspond to larger magnitudes of the cumulants. η 2 < 0 corresponds to the first order phase transition region and η 2 > 0 corresponds to the crossover region. η 1 = 0 refers to the phase transition line. In the subfigure κ 4 , the dotted blue lines mark the boundary where κ 4 changes the sign; while as a comparison, the dashed purple lines are the boundary in the thermodynamic limit case. The bottom-right subfigure shows the setup of the dynamical evolution in Sec. 5, and visualize the memory effects with different relaxation rates. The grey lines are the contours of µ and T in the η 1 -η 2 plane. The dynamical evolution are set to be along fixed chemical potentials, while the blue and red solid line are the starting and hypothetical freeze-out line for the dynamical evolution. The evolution direction is marked by arrows. The dashed red (green, blue) line, on which the equilibrium κ 1 equals to the dynamical κ 1 on the solid red line (hypothetical freeze-out line) with relaxation rate τ rel = 0.05(0.1, 0.2 fm), represents the effective freeze-out line due to the memory effects of dynamics.
ergy. The Fokker-Plank equation for the dynamical probability distribution (denoted by P[σ; t]) of the σ field is [31] where Ω[σ; t] ≡ − T V lnP[σ; t] is the dynamical free energy density and Ω 0 [σ] is the equilibrium free energy density at the point (µ(t), T (t)) (see Eq. (1)). The parameter τ eff is the effective relaxation rate. m σ is the equilibrium mass of the σ field, which is defined as where σ 0 is the global minimum of Ω 0 [σ]. The following calculation assumes that the dependence of the relaxation rate τ eff on the equilibrium correlation length ξ eq (µ, T ) satisfies where ξ eq = m −1 σ [31]. Here τ rel and ξ ini are the initial relaxation rate and the initial equilibrium correlation length, respectively. The value z = 3 is given by the dynamical critical exponent of Model H [68,69].
The time evolution equation of Ω[σ; t] is deduced from Eq. (16), and we obtain By denoting Ω[σ; t] = i α i (t)σ i /i, the time evolution of the coefficient α i (t) becomes, where α 0 j ( j = 1, · · · , 4) are the coefficients for the equilibrium free energy density Ω 0 [σ]. Note that different coefficients of the σ field couple with each other in the time-evolution equations. As a result, the terms of σ i with power i > 4 emerge automatically after evolution, even though they vanish at the initial setting. As before, by assuming the contributions from higher power terms are negligible, the free energy is cut off till the fourth order. By setting up the parameters of the dynamical system, we can numerically solve the coupled equations in above, and study the dynamical behaviors of the physical quantities like the cumulants.

non-equilibrium cumulants
As we have obtained the time evolution equation (20) of the QCD free energy, we can study the corresponding dynamical cumulants. For simplicity, in the following calculations, we suppose that both the chemical potential and the volume (V = 10 3 fm 3 ) are fixed during the time evolution 1 . The temperature is assumed to decrease as a function of time:  In Fig. 4, we plot the evolution of the cumulants for different relaxation rates, starting from an initial temperature T ini = 185 MeV. The monotonicity of the equilibrium cumulants are exactly memorized by the non-equilibrium cumulants, and reproduced in the later evolution process. As the increase of relaxation rate, the peaks and dips of the non-equilibrium cumulants are suppressed, and the typical structure for each order of the cumulants appears in a later time. Note that in the current finite size system, the memory effects shown in Fig. 4 is similar for different chemical potentials for both the crossover and the first order phase transition cases.
Next, we study the dynamical result of the σ's cumulants at a hypothetical freeze-out line. The initial temperature and the hypothetical freeze-out temperature are supposed to be about 5.05 MeV above and 3.37 MeV below the phase transition temperature, respectively (i.e. η 1 = 1.5 × 10 5 MeV 3 for the initial state and η 1 = −1.0 × 10 5 MeV 3 for the final state, marked as blue and red lines in the bottom-right subfigure in Fig. 3). In Fig. 5 we present the non-equilibrium cumulants with respect to the chemical potential. With different relaxation rates, the non-equilibrium cumulants at large chemical potentials (mainly on the first order side) significantly deviate from the equilibrium cumulants, and develop rich structures. Especially, for the blue line with τ rel = 0.2 fm, κ 3 change its sign, and κ 4 presents non-monotonic behaviors compared to the equilibrium ones. These rich structures are suggestive for the understanding of the STAR data [27].
The dynamical behaviors of the cumulants can be understood from the equilibrium cumulants by considering the memory effects. Let us define an effective freeze-out temperature (for each chemical potential) where the equilibrium expectation of the σ field, κ 1 , equals to the dynamical κ 1 on the freezeout point. The effective freeze-out lines for different relaxation rates are shown in the bottom-right subfigure of Fig. 3 (The dashed red, green and blue lines). From the equilibrium cumulants κ 2 -κ 4 on the effective freeze-out lines, we can approximately read the sign of the non-equilibrium cumulants in Fig. 5. Specifically, for κ 3 , when the relaxation rate is large, say τ rel = 0.2 fm, the effective freeze-out line in the crossover region is basically at the left side of the phase transition line (η 1 < 0), while in the first order phase transition region, it is still at the right side of the phase transition line (η 1 > 0), leading to a positive dynamical κ 3 (see κ 3 in Fig. 3). For κ 4 , when τ rel = 0.1 fm, the effective freeze-out line falls into the negative κ 4 region; when τ rel = 0.2 fm (τ rel = 0.05 fm), the effective freeze-out line crosses the right sign-change line of the κ 4 in η 1 -η 2 plane (see κ 4 in Fig. 3). In general, the memory effects in the first order phase transition region is more significant than that in the crossover region (i.e. the cumulants record earlier information at larger µ), which is because in the first order phase transition region, the barrier in the free energy strongly delays the dynamical evolution of the σ field.

Discussion
In summary, we parameterize the Landau free energy of the χPT by mapping it to the Ising free energy. With the parametric free energy, we study the equilibrium cumulants in a finite size system. The volume of the system significantly influences the fluctuations of the σ field. In the equilibrium case, we find that for a typical QGP fireball size with volume V = 10 3 fm 3 , the probability distribution of the σ field is broadened near the phase transition line, leading to enhancements of the fluctuations and rounding of the discontinuities on the first order phase transition side. Moreover, the fourth order cumulant κ 4 of the σ field is universally negative in the phase transition region on both the crossover and the first order phase transition side. Compared to the crossover region, all cumulants in the first order phase transition region are enhanced, due to the two-peak shape of σ's distribution in the finite system (where the barrier of the free energy density is of the order of T/V). Utilizing the Fokker-Plank equation, we derive the real-time evolution of the parametric free energy, and further, by setting the cooling down process of the system, we analyze the time evolution of non-equilibrium cumulants at different phase transition scenarios and non-equilibrium cumulants on the hypothetical freezeout line. We find that earlier information about the cumulants are recorded in the first order phase transition region, compared to the crossover region. Like the earlier studies, the dynamical cumulants can be different from the equilibrium ones from both the magnitude and the sign, and can be nonmonotonic on the hypothetical freeze-out line.
Note that in our study, we adopt the zero mode approximation which assumes the order parameter is uniform in space. The nonzero modes can be further taken into account if one employs the Ginzburg-Landau free energy. Moreover, due to the flexibility of our parameterization, the parametric free energy and its dynamical evolution can be easily combined with a realistic dynamical model (like that in chiral hydrodynamics [38,39,40,41] or Hydro+ [42]), to study the full dynamical process with phase transition. Now compare our method with the former studies on the non-equilibrium cumulants in Refs. [31,32,33,34]. In Ref. [31], they presented pioneering studies on the real-time evolution of non-Gaussian cumulants derived from the Fokker-Plank equation. The QCD equation of state near the critical point in the thermodynamic limit is inspired from the Ising model and the parameterization is realized by linearly mapping QCD variables (µ, T ) to the Ising variables (r, h). The resulted equation of state and the corresponding equilibrium cumulants could contain the quantum correction from the fluctuations of the order parameter. However, the volume effects on the equation of state and cumulants are fully omitted, which is valid when the correlation length is significantly smaller than the size of the system. In this case, the discontinuities of the equilibrium distribution function of the σ field hinder the simulation of critical dynamics with the Fokker-Plank equation in the first order phase transition region, which assumes the slight deviations of dynamical probability distribution from the equilibrium one. As shown in our study, when the size is small enough, the finite size effects becomes significant especially in the first order phase transition region. The rounding effects justifies the application of the Fokker-Plank equation in the first order phase transition region. In Refs. [32,33,34], the dynamical evolution of cumulants on different phase transition scenarios were studied based on event-by-event simulations of Langevin equation. The effective potential of σ is obtained by evaluating the linear sigma model. The evolution takes account of the fluctuation effects (long wavelength modes) in the real space. However, the shortcoming of this approach is the location of the critical point is fixed by the model parameters of the linear sigma model, which is even far below the chemical freeze-out line determined by the statistical model and experimental data. Thus the dynamical equations of the critical fluctuations in Ref. [32,33] could not be combined with the hydrodynamic equations, and describe the realistic experimental process. In recent studies, the parameterization of QCD equation of state with nonlinear corrections (consistent with the lattice results) have been extended from the critical region to a broader region in the µ-T plane [50,51,52]. Similar further extension of the parameterization of the free energy density will be important for studying the dynamic evolution of critical fluctuations at RHIC.