# Ising parameterization of QCD Landau free energy and its dynamics

###### Abstract

We present a general linear parameterization scheme for the QCD Landau free energy in the vicinity of the critical point of chiral phase transition in the - plane. Based on the parametric free energy, we show that, due to the finite size effects, the regions of fluctuations of the order parameter (i.e. the field) are broadened, and the discontinuities of the first order phase transition are smoothed. Meanwhile, the kurtosis of the field is universally negative around the critical point. Using the Fokker-Plank equation, we derive the dynamical corrections to the free energy. The dynamical cumulants of the field on the freeze-out line, record earlier information in the first order phase transition region as compared to the crossover region. The typical behavior of dynamical cumulants can be understood from the equilibrium cumulants by considering the memory effects.

## I Introduction

The QCD phase structure has been of intensive theoretical and experimental interest for decades Cabibbo:1975 (); Baym:2002 (); Fukushima:2011 (); Aggarwal:2010cw (). In various extreme conditions, the QCD matter could exist as quark-gluon plasma (QGP), hadronic resonance gas (HRG), or color superconductor, etc. A conjecture for the phase diagram in the - plane is summarized in Nuclear:2008 (). Of great interest is the chiral phase transition (PT) between the hot QGP phase and the HRG phase stocker:1986 (); TDLee:1974 (); Hatsuda:1985 (). For small baryon chemical potential, , lattice simulations show that the transition is a broad crossover Brown:1990 (); Aoki:2006we (); Aoki:2009sc (); Karsch:2003jg (). For large , the sign problem of lattice QCD prevents full ab initio simulations. Hence, the information about the order of the phase transition can only be obtained through effective models, like the Nambu-Jona-Lasinio model Klevansky:1992qe (); Fu:2007xc () and quark-meson model Schaefer:2007 (), or non-perturbative methods like Functional Renormalization Group Berges:2000ew () and Dyson-Schwinger equation Roberts:1994dr (). These models show that the transition is of first order for large . These phenomenological analysis predict the existence of a chiral critical point.

In the vicinity of the critical point, the correlation length of the field is largely increased, leading to intense fluctuations of the net proton production according to perturbation theory Stephanov:2008qz (). The related observables are presently measured by the STAR collaboration at the Relativistic Heavy Ion Collider (RHIC) Adamczyk:2013dal (); Luo:2015ewa (). However, these current experimental data can not even qualitatively be explained by theoretical calculations Jiang:2015hri (); Mukherjee:2015swa (); Bluhm:2016byc (), due to the complexities of a realistic dynamical modeling of the experimental situation Paech:2003fe (); Nahrgang:2011mg (); Nahrgang:2011vn (); Herold:2013bi (); Herold:2013qda (); Stephanov:2017ghc (). A main component of the dynamical modeling is the parameterization of the QCD free energy or the equation of state near the QCD critical point.

As the critical behaviors of the field fall into the same universality class as the 3D Ising model Gavin:1994 (), one can obtain the parametric equation of state by mapping the QCD parameters onto the Ising variables Guida:1997 (); Berdnikov:1999ph (); Nonaka:2005 (); Stephanov:2011 (). Combined with the dynamical evolution function, there are various interesting dynamical effects presented, like memory effects Mukherjee:2015swa () and universal off-equilibrium scaling behaviors Berdnikov:1999ph (); Mukherjee:2016kyu (). However, the exact mapping from to is still unknown, which leads to large uncertainties in the static and dynamical study. Also, the parametric equation of state is obtained in the thermodynamic limit Cardy:1996 (); Kos:2016 (); Komargodski:2017 (), which may not be applicable to the RHIC physics with a finite size Yang:1952 (); Spieles:1997ab ().

In the present paper, we develop a general linear parametric scheme for the QCD Landau free energy in the vicinity of the critical point, by comparing with the free energy of the Ising model. The method can include finite size effects naturally. Using the parametric free energy, we study the finite size effects and the equilibrium cumulants of the field. The dynamical corrections to the free energy and the cumulants are evaluated by employing the Fokker-Plank equation.

Let us point out that for a finite size system, the discontinuities of the first order phase transition (FOPT) are smoothed, and the distribution of the field is broadened. Consequently, the kurtosis of the field is universally negative, both on the crossover side and on the FOPT side. For the dynamical evolution, the cumulants record earlier information in the FOPT region, as compared to the crossover region.

The article is structured as follows. In section II, we develop the Ising mapping of the QCD Landau free energy. The linear parameterization of the free energy is realized. In section III, the dynamical free energy is obtained through the Fokker-Plank equation. The static finite size effects are discussed in section IV and the dynamical cumulants are evaluated in section V. Finally, we give more discussions.

## Ii Ising parameterization

According to the Landau theory of phase transitions, the free energy in the critical region is supposed to be analytic and obey the symmetry of the Hamiltonian. The Landau free energy density of the PT Pisarski:1984 (); Gell-Mann:1960 () can be written in terms of the field as

(1) |

with the Taylor expansion of the done up to the fourth order. The constant term is omitted here because it is irrelevant to the structure of the QCD phase diagram. A zero-momentum mode approximation of the field is assumed. The distribution of the critical fluctuations is described by the probability function Stephanov:2008qz (), where is the volume of the system. In the chiral limit, Petropoulos:1999 (), while for the physical world, a finite term, , is introduced to handle the explicit chiral symmetry breaking of the quark masses. The cubic term, , emerges only after the renormalization contributed from the high-momentum modes of the field. The coefficient of the quartic term, , is supposed to be positive, in order to sustain the stability of the system.

In the thermodynamic limit (), the phase structure is fully determined by the global minimum of the free energy (1). The information on the PT is presented by performing a translation transformation for the sigma field, , with

(2) |

Consequently, the cubic term in eq. (1) can be eliminated, and we obtain the translated free energy density

(3) |

which has exactly the same form as the free energy density of the Ising model. The redefined coefficients, , can be expressed as functions of the original coefficients and the shift by

(4) | |||||

(5) | |||||

(6) |

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 - plane into four parts, according to the sign of and . As shown in Fig. 1, on the phase transition line, , the sign of determines the order of the phase transition. In regions II and IV, , two minima coexist in the coexistence region, which is determined by . The physical vacuum is the lower one, which is the typical character of the FOPT. On the other part of the - plane, , only one minimum exists, and the system undergoes a continuous phase transition (crossover) as crosses its zero point. The sign of determines the location of the global minimum. In general, in the upper regions I and II, with , the global minimum is at (i.e. ); while in the lower regions III and IV, with the global minimum is at (i.e. ). At the critical point , both and vanish.

To parameterize the coefficients in linearly in the - plane, it is convenient to define two unit vectors and . The vector is parallel to the tangent of the phase transition curve at the critical point. Note that and are not necessarily orthogonal. The angle between them should be determined by experiment. Around the critical point, by projecting the vector onto the perpendicular vector of and that of , respectively, the coefficients are linearly parameterized as:

(7) | |||||

(8) | |||||

(9) |

With these projections, the sign of and in the different regions (I-IV) can be correctly expressed by constraining all the constants to (). Along the phase transition curve, , the correlation length is , which diverges at the critical point, and is finite on both the FOPT side and the crossover side Fukushima:2011 (). In the critical region, the variable is parameterized as a constant, , with the zeroth order approximation. Then, the original free energy (eq.(1)) can also be linearly parameterized, through the relations (2) and (4)-(6).

In the current Ising parameterization of the QCD free energy, the coefficients and in the translated free energy (3) correspond exactly to the Ising variables and , respectively. The parameterization used here is universal. Hence, it can also be applied to the (hadronic) gas-liquid phase transition Vovchenko:2015pya (); Vovchenko:2016rkn () or to other similar scenarios, by replacing the field with the corresponding order parameter. Moreover, it can be generalized to finite size systems (see the section ‘finite size effects’).

Through the transformation of eq. (1) to eq. (3), we present the connection between the QCD free energy and the Ising free energy. Different from the former approaches on the QCD equation of state Guida:1997 (); Berdnikov:1999ph (); Nonaka:2005 (); Stephanov:2011 (), we provide a new parametric approach on the QCD free energy. Within linear approximation, it makes the direct mapping from the curved QCD phase transition lines in - plane to the Ising phase transition lines in the - plane.

Through the transformation of eq. (1) to eq. (3), we present the similarity between the QCD free energy and the Ising free energy. It provides a new mapping method different from the former ones based on the QCD equation of state Guida:1997 (); Berdnikov:1999ph (); Nonaka:2005 (); Stephanov:2011 (). The new mapping is rather direct. It makes the mapping from the curved QCD phase transition lines in - plane to the Ising phase transition in the - plane within linear approximation.

## Iii Dynamical free energy

With the parameterized equilibrium free energy, we now develop the dynamical equations. The Fokker-Plank equation for the dynamical probability distribution (denoted by ) of the field is Mukherjee:2015swa ()

(10) |

Here, is the dynamical free energy density. is the equilibrium free energy density at (see eq. (1)). The parameter is the effective relaxation rate. is the equilibrium mass of the field, defined as , where is the global minimum of . The following calculation assumes that, the dependence of the relaxation rate on the equilibrium correlation length satisfies Mukherjee:2015swa (). Here and are the initial relaxation rate and the initial equilibrium correlation length, respectively. The value is given by the dynamical critical exponent of Model H Hohenberg:1977 (); Son:2004iv ().

The evolution equation of is obtained from eq. (10),

(11) | |||||

Denoting , the evolution of each coefficient becomes,

(12) | |||||

Here is the coefficient of . The terms, , with power emerge automatically during the time evolution from the above coupled equations. As before, the free energy is cut at the order , assuming that the contributions from higher power terms are negligible.

## Iv Finite size effects

The discussion of the free energy density (1) is based on the assumption of the thermodynamic limit. In reality, the size of the heavy ion system is finite. This has a direct influence on the renormalization process. Consequently, the coefficients in eq. (1) depend on the volume brezin1985npb (). Hereafter, we will not focus in detail on these finite size corrections to the free energy density. Rather, we point out that the parameters , , , , , , , and are volume dependent. The numerical calculations, are done for the following parameters: , , , , , , and . Here, is set for simplicity. The values of these parameters are estimated so that the cumulants are .

Even without considering its effects on the free energy density, the finite volume can still influence the fluctuations of the field through the probability function . Fig. 2 shows the cumulants ^{1}^{1}1The cumulants of the field are defined as , , , and , with the moments . of the field with respect to the temperature for different volumes. For a large enough volume (say ), the width of the field’s distribution is suppressed, and the field locates only at the global minimum. While for a smaller system, the equivalent ‘temperature’ is increased, which leads to a larger variance for the distribution of the field. The expectation value of the order parameter, , deviates from the global minimum of the free energy. It becomes smooth even for the FOPT, as the field occupies the two minima simultaneously around the phase transition temperature. Similar finite-size rounding effect of a FOPT can also be found in Ref. Imry1980prb (); Spieles:1997ab (); Zabrodin:1998vt (). For a small system, the perturbation theory around the global minimum Stephanov:2008qz (); Jiang:2017mji (); Jiang:2017fas () does not work well anymore.

Fig. 3 presents the contours of the cumulants in the - plane for , which confirm again large cumulants (-) in the FOPT region. Similar to Stephanov:2008qz (), changes its sign after crossing the phase transition line. On the other hand, the kurtosis () is generally negative around the critical point. This is due to the fact that, both the two-peak shape and the intermediate shape between Gaussian and two-peak shape are broader than the Gaussian Stephanov:2011 (), as long as the system is near the phase transition temperature. This negative kurtosis is in striking contrast to the results in the thermodynamic limit. There, the kurtosis is negative only in the crossover region (the line of the change of the sign is marked by the dashed purple line in subfigure ), as only the global minimum is occupied in the FOPT region Stephanov:2011 ().

## V Dynamical cumulants

The evolution equation (12) is employed to study the dynamical cumulants. Along the trajectory, both the chemical potential and the volume are supposed to be fixed disc (). The evolution of the temperature satisfies , where is the initial temperature, is the reference time and Mukherjee:2015swa (). In Fig. 4, we plot the evolution of the cumulants for different relaxation rates, starting from the initial temperature . The memory effect shown in Fig. 4 is similar for all chemical potentials around . The peaks and dips of the cumulants are suppressed during the evolution with increasing relaxation rate. Fig. 5 shows the dynamical cumulants with respect to the chemical potential. The initial and the freeze-out temperature are supposed to be about above and below the phase transition temperature, respectively (i.e. for the initial state and for the final state, see the lower r.h.s. subfigure in Fig. 3). For different relaxation rates, the dynamical cumulants on the FOPT side significantly deviate from the equilibrium cumulants, and present rich structures. Especially, for , the change of sign for and the presence of non-monotonic behavior for compared to equilibrium ones are suggestive for the understanding of the STAR data Luo:2015ewa ().

The dynamical behavior of the cumulants can be understood from the equilibrium cumulants by considering the memory effects. Let us define an effective temperature where the equilibrium expectation of the sigma field, , equals to the dynamical on the freeze-out line. The effective temperature lines for different relaxation rates are shown in the lower r.h.s. subfigure of Fig. 3. The equilibrium cumulants - on the effective temperature lines, approximately reflect the sign of the dynamical cumulants in Fig. 5. Specifically, for , when the relaxation rate is large, say , the effective temperature in the crossover region is basically below the phase transition line (), while in the FOPT region, it is still above the phase transition line (), leading to a positive dynamical . For , when , the effective temperature falls into the negative kurtosis region; when (), the effective temperature line crosses the right (left) sign-change line of the kurtosis. In general, the memory effect in the FOPT region is more significant than that in the crossover region (i.e. the cumulants record earlier information). This is due to the barrier in the free energy in the FOPT region.

## Vi Discussion

Now we compare our method with the former studies on the dynamical cumulants in Mukherjee:2015swa (); Jiang:2017mji (); Jiang:2017fas (). In Mukherjee:2015swa (), the dynamical evolution of cumulants are derived from the Fokker-Plank equation. The parametric QCD equation of state near the QCD critical point in the thermodynamic limit is taken as input. However, the mapping from Ising variables to the QCD variables is not clear, as pointed out in their paper. The advantage of this method is that it contains the quantum fluctuation correction to the critical indexes. However, the method can not be extended directly to the FOPT region when the finite size effect becomes significant. In Jiang:2017mji (); Jiang:2017fas (), the dynamical evolution of high order cumulants are studied based on event-by-event simulations of Langevin equation. The effective potential of is obtained by evaluating the linear sigma model. It takes account of the fluctuation effect (long wavelength modes) in the real space. However, the critical point in the model depends on the number of flavors of the quarks and the coupling strength of the interaction, and locates in the hadron resonance gas region, below the statistical freeze-out line.

In our calculations, the Ising mapping is built between the QCD free energy and Ising free energy by using their similarities. The dynamical free energy is derived from the Fokker-Plank equation. The corresponding dynamical cumulants are calculated based on the dynamical free energy. In this paper, we use the zero mode approximation, which is intrinsic a mean field approximation. The nonzero modes can be further taken into account if we employ the Ginzburg-Landau free energy. Thus the mapping is compatible with Langevin equation. This method can be easily extended to finite size system, and is also applicable in the FOPT region. Furthermore, the mapping is model independent and the parameterization of free energy is universal within the 3D Ising model universality class.

In summary, we have provided a new parametric approach to the QCD free energy near the critical point by comparing it to the Ising model. The parameterization can be applied in all phase transitions within the same universality class. The volume of the system significantly influences the fluctuations of the field. A small volume smoothes the discontinuities in the FOPT region. In addition, compared to the crossover region, the cumulants in the first order phase transition region are enhanced, due to the broadening of ’s distribution. The kurtosis is generally negative around the critical point. We also study the dynamical evolution of the free energy and the related cumulants. We find that earlier information about the cumulants are recorded in the FOPT region, rather than in the crossover region. The sign of the dynamical cumulants are reflected in general by the equilibrium cumulants, after considering the memory effects.

Note that the current parameterization of the QCD free energy is linear, it works in the region close to the critical point. In a wider region, higher order expansions of (, T) should be taken into account. The parametric free energy and its dynamical evolution can be easily employed in a realistic dynamical model (like chiral hydrodynamicsNahrgang:2011mg (); Nahrgang:2011vn (); Herold:2013bi (); Herold:2013qda () or Hydro+Stephanov:2017ghc ()) in RHIC, with tunable location of the critical point and controllable range of the critical region as input. The realistic dynamical modeling with this parametric QCD free energy will be done in future studies.

## Vii Acknowledgments

L. Jiang thanks Volodymyr Vovchenko, Shanjin Wu and Huichao Song for reading the manuscript and comments. H. Stöcker acknowledges the support through the Judah M. Eisenberg Laureatus Chair at Goethe University. This work is supported by the GSI Helmholtzzentrum für Schwerionenforschung and by the Helmholtz International Center for the Facility for Antiproton and Ion Research (HIC for FAIR) within the framework of the Landes-Offensive zur Entwicklung Wissenschaftlich-Oekonomischer Exzellenz (LOEWE) program launched by the State of Hessen.

L. Jiang and J.-H. Zheng contributed equally to this work.

## References

- (1) N. Cabibbo and G. Parisi, Phys. Lett. B 59, 67 (1975).
- (2) G. Baym, Nucl. Phys. A 956, 1 (2016).
- (3) K. Fukushima and T. Hatsuda, Rep. Prog. Phys. 74, 014001 (2011).
- (4) M. M. Aggarwal et al. [STAR Collaboration], arXiv:1007.2613.
- (5) The DOE/NSF Nuclear Science Advisory Committee, arXiv:0809.3137 (2008).
- (6) H. Stöcker and W. Greiner, Phys. Reps. 137, 277 (1986).
- (7) T. Hatsuda and T. Kunihiro, Phys. Rev. Lett. 55, 158 (1985).
- (8) T. D. Lee, and G. C. Wick, Phys. Rev. D, 9, 2291 (1974).
- (9) F. R. Brown, et al, Phys. Rev. Lett. 65, 2491 (1990).
- (10) Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz and K. K. Szabo, Nature 443, 675 (2006).
- (11) Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz, S. Krieg and K. K. Szabo, JHEP 0906, 088 (2009).
- (12) F. Karsch and E. Laermann, In *Hwa, R.C. (ed.) et al.: Quark gluon plasma* 1-59 [hep-lat/0305025].
- (13) S. P. Klevansky, Rev. Mod. Phys. 64, 649 (1992); K. Fukushima, Phys. Lett. B 591, 277 (2004).
- (14) W.J. Fu, Z. Zhang and Y. X. Liu, Phys. Rev. D 77, 014006 (2008); L.J. Jiang, X.Y. Xin, K.L. Wang, S.X. Qin and Y.X. Liu, Phys. Rev. D 88, 016008 (2013).
- (15) B.-J. Schaefer and J.Wambach, Nucl. Phys. A 757, 479 (2005); B.-J. Schaefer, J. M. Pawlowski, and J. Wambach Phys. Rev. D 76, 074023 (2007); T. K. Herbst, J. M. Pawlowski and B.-J. Schaefer, Phys. Lett. B 696, 58 (2011).
- (16) J. Berges, N. Tetradis and C. Wetterich, Phys. Rept. 363, 223 (2002).
- (17) C. D. Roberts and A. G. Williams, Prog. Part. Nucl. Phys. 33, 477 (1994); S.X. Qin, L. Chang, H. Chen, Y.X. Liu and C.D. Roberts, Phys. Rev. Lett. 106, 172301 (2011).
- (18) M. A. Stephanov, Phys. Rev. Lett. 102, 032301 (2009).
- (19) L. Adamczyk et al. [STAR Collaboration], Phys. Rev. Lett. 112, 032302 (2014).
- (20) X. Luo [STAR Collaboration], PoS CPOD 2014, 019 (2014).
- (21) L. Jiang, P. Li and H. Song, Phys. Rev. C 94, 024918 (2016).
- (22) S. Mukherjee, R. Venugopalan and Y. Yin, Phys. Rev. C 92, 034912 (2015).
- (23) M. Bluhm, M. Nahrgang, S. A. Bass and T. Schaefer, Eur. Phys. J. C 77, 210 (2017).
- (24) K. Paech, H. Stöcker and A. Dumitru, Phys. Rev. C 68, 044907 (2003).
- (25) M. Nahrgang, S. Leupold, C. Herold and M. Bleicher, Phys. Rev. C 84, 024912 (2011).
- (26) M. Nahrgang, C. Herold, S. Leupold, I. Mishustin and M. Bleicher, J. Phys. G 40, 055108 (2013).
- (27) C. Herold, M. Nahrgang, I. Mishustin and M. Bleicher, Phys. Rev. C 87, 1, 014907 (2013).
- (28) C. Herold, M. Nahrgang, I. Mishustin and M. Bleicher, Nucl. Phys. A 925, 14 (2014).
- (29) M. Stephanov and Y. Yin, arXiv:1712.10305 [nucl-th].
- (30) S. Gavin, A. Gocksch and R. D. Pisarski, Phys. Rev. D 49 3079 (1994).
- (31) R. Guida and J. Zinn-Justin, Nucl. Phys. B 489, 626 (1997).
- (32) B. Berdnikov and K. Rajagopal, Phys. Rev. D 61, 105017 (2000).
- (33) C. Nonaka and M. Asakawa, Phys. Rev. C 71 044904 (2005).
- (34) M. A. Stephanov, Phys. Rev. Lett. 107, 052301 (2011).
- (35) S. Mukherjee, R. Venugopalan and Y. Yin, Phys. Rev. Lett. 117, 222301 (2016).
- (36) J. Cardy, Scaling and Renormalization in Statistical Physics, Cambridge University Press, ISBN 978-0-521-49959-0 (1996).
- (37) F. Kos, D. Poland, D. Simmons-Duffin and A. Vichi, JHEP 2016: 36 (2016).
- (38) Z. Komargodski, and D. Simmons-Duffin, J. Phys. A: Mathematical and Theoretical, 50 154001 (2017).
- (39) C. N. Yang and T. D. Lee, Phys. Rev. 87, 404 (1952); T. D. Lee and C. N. Yang, Phys. Rev. 87, 410 (1952).
- (40) C. Spieles, H. Stöcker and C. Greiner, Phys. Rev. C 57, 908 (1998).
- (41) R. D. Pisarski and F. Wilczek, Phys. Rev. D 29 338 (1984).
- (42) M. Gell-Mann and M. Levy, Nuovo Cimento 16 705 (1960).
- (43) N. Petropoulos, J. Phys. G: Nuclear and Particle Physics 25 2225 (1999).
- (44) V. Vovchenko, D. V. Anchishkin, M. I. Gorenstein and R. V. Poberezhnyuk, Phys. Rev. C 92, 054901 (2015).
- (45) V. Vovchenko, M. I. Gorenstein and H. Stoecker, Phys. Rev. Lett. 118, 182301 (2017).
- (46) P. C. Hohenberg, and B. I. Halperin, Rev. Mod. Phys. 49, 435 (1977).
- (47) D. T. Son and M. A. Stephanov, Phys. Rev. D 70, 056001 (2004).
- (48) E. Brezin and J. Zinn-Justin, Nucl. Phys. B 257, 867 (1985).
- (49) Y. Imry, Phys. Rev. B 21, 2042 (1980).
- (50) E. E. Zabrodin, L.V. Bravina, H. Stöcker and W. Greiner, Phys. Rev. C 59, 894 (1999).
- (51) L. Jiang, S. Wu and H. Song, Nucl. Phys. A 967, 441 (2017).
- (52) L. Jiang, S. Wu and H. Song, EPJ Web Conf. 17, 11600 (2018).
- (53) The time of evolution around the phase transition curve turns out to be rather brief in this numerical calculation. Hence, the change of the volume during the expansion is omitted here.