Baryon number fluctuations in chiral effective models and their phenomenological implications
BBaryon number fluctuations in chiral effective models and their phenomenologicalimplications
G´abor Andr´as Alm´asi, ∗ Bengt Friman, and Krzysztof Redlich
2, 3 GSI Helmholtzzentrum f¨ur Schwerionenforschung, D-64291 Darmstadt, Germany University of Wroc(cid:32)law, Institute of Theoretical Physics, PL 50-204 Wroc(cid:32)law, Poland ExtreMe Matter Institute EMMI, GSI, D-64291 Darmstadt, Germany
We study the critical properties of net-baryon-number fluctuations at the chiral restoration tran-sition in a medium at finite temperature and net baryon density. The chiral dynamics of quantumchromodynamics (QCD) is modeled by the Polykov-loop extended Quark-Meson Lagrangian, thatincludes the coupling of quarks to vector meson and temporal gauge fields. The Functional Renor-malization Group is employed to properly account for the O (4) criticality at the phase boundary.We focus on the properties and systematics of ratios of the net-baryon-number cumulants χ nB , for1 ≤ n ≤
6, near the phase boundary. The results are presented in the context of the recent ex-perimental data of the STAR Collaboration on fluctuations of the net proton number in heavy ioncollisions at RHIC. We show that the model results for the energy dependence of the cumulant ratiosare in good overall agreement with the data, with one exception. At center-of-mass energies below19 . O (4) and Z (2) criticality. We assess the influence of modelassumptions and in particular of repulsive vector-interactions, which are used to modify the locationof the critical endpoint in the model, on the cumulants ratios. Finally, we discuss a possibility totest to what extent the fluctuations are affected by nonequilibrium dynamics by comparing certainratios of cumulants. I. INTRODUCTION
It is well established within lattice QCD (LQCD), thatat vanishing and small chemical potential, µ B , stronglyinteracting matter undergoes a smooth crossover fromthe hadronic phase to a quark-gluon plasma [1–4]. More-over, arguments are given, that in the limit of massless u and d quarks, the crossover transition becomes a gen-uine second-order chiral phase transition, belonging tothe O (4) universality class [5, 6]. The nature of this tran-sition at higher net baryon densities is still not settled byfirst principle LQCD studies, owing to the sign problem.However, in effective models of QCD, it is found that atsufficiently large µ B , the system can exhibit a first orderchiral phase transition. The endpoint of this conjecturedtransition line in the ( T, µ B )- plane, is the chiral criticalendpoint (CEP) [7–11]. At the CEP, the system exhibitsthe 2 nd order phase transition, which belongs to the Z (2)universality class [12].Various methods have been developed to circumventthe sign problem and perform LQCD calculations at non-vanishing µ B [13–15]. However, so far these methods arerestricted to rather small values of chemical potential, µ B (cid:28) T c , where T c = 154(9)MeV is the chiral crossovertemperature at vanishing net baryon density [3, 4].Due to the restriction of present LQCD calculations tosmall net baryon densities, effective models that belongto the same universality class as QCD, e.g., the Polyakov-loop extended Nambu-Jona-Lasinio [16–24] and Quark-Meson (PQM) models [25–33], have been employed to ∗ [email protected] study the chiral phase transition for a broad range ofthermal parameters. Moreover, there are ongoing stud-ies using functional methods that focus on systematicimprovements of such effective models, to achieve a morequantitative matching with LQCD, and eventually toperform reliable calcuations of the properties of QCDmatter at nonzero net baryon density [34–36].One of the strategic goals of current experimental andtheoretical studies of chiral symmetry restoration in QCDis to unravel the phase diagram of strongly interactingmatter and to clarify whether a chiral CEP exists. A ded-icated research program at RHIC, the beam energy scan(BES), has been established to explore these issues in col-lisions of heavy ions at relativistic energies [37]. By vary-ing the beam energy at RHIC, the properties of stronglyinteracting matter in a broad range of net baryon den-sities, correponding to a wide range in baryon chemi-cal potential, 20 MeV < µ B <
500 MeV [38–41], can bestudied. In particular, fluctuations of conserved chargesare considered as relevant probes of the phase dia-gram [23, 24, 42–51]. These are experimentally acces-sible and reflect the criticality of the chiral transition.Fluctuations of the net baryon number are particularyinteresting, owing to a direct connection to critical scal-ing near the chiral phase boundary [46, 51].First data on net-proton-number fluctuations, whichare used as a proxy for fluctuations of the net baryonnumber, in heavy-ion collisions have been obtained bythe STAR Collaboration at RHIC energies [52–54]. TheSTAR data on the variance, skewness and kurtosis ofnet proton number, are intriguing and have stimulatedlively discussions on their physics origin and interpreta-tion. Recently, first results on the mean and variance ofthe net-proton distribution were obtained in heavy ion a r X i v : . [ h e p - ph ] M a r collisions at the LHC energy by the ALICE Collabora-tion [55].In the following, we focus on the properties and sys-tematics of the cumulants of net-baryon-number fluctu-ations near the chiral phase boundary. Our studies areperformed based on the Polykov-loop extended Quark-Meson Lagrangian which includes couplings of quarks tovector and temporal gauge fields. To account for the O (4) and Z (2) critical fluctuations near the phase bound-ary, we employ the Functional Renormalisation Group(FRG) [56–59]. We present contour plots of ratios ofcumulants, involving the mean, the variance, skewnessand kurtosis in the ( T, µ B )- plane. These ratios obtainedon the phase boundary and on a freeze-out line deter-mined by fitting the skewness ratio, following Ref. [60],are confronted with the corresponding experimental val-ues of the STAR Collaboration. We explore the influenceof the CEP and the repulsive interactions on the fluctu-ation observables. The calculations are performed withtwo different initial conditions and two Polyakov-loop po-tentials, to assess the model dependence of the results.We find that the energy dependence of ratios of low-order cumulants, χ nB with n = 1 , ,
3, are in good agree-ment with the data, whereas for the ratio involving thekurtosis, the model results differ from the data at en-ergies below √ s (cid:39)
20 GeV, corresponding roughly tothe top SPS energy. At these low energies, the latterratio increases strongly, while the model results, whichembody O (4) as well as Z (2) criticality, differ substan-tially from the data. Finally, we discuss possible caveats,which could undermine our conclusions and assess theuncertainties of the model.The paper is organized as follows. In Sec. II we formu-late the model and the FRG method, which is employedto compute cumulants of net baryon number. In Sec. IIIwe present results on the characteristics of susceptibilityratios as functions of temperature and chemical potentialand study, along the phase boundary, the dependence ofthese ratios on the vector interaction. Moreover, in thissection we confront the model results with the STARdata on net-proton-number fluctuations. Sec. IV con-tains summary and conclusions. II. THE POLYAKOV–QUARK-MESON MODEL
The PQM model is a low energy effective approxi-mation to QCD formulated in terms of the light quark q = ( u, d ) as well as scalar and the pseudoscalar meson φ = ( σ, (cid:126)π ) fields. The quarks are coupled to the back-ground Euclidean gluon field A µ , with vanishing spatialcomponents, which is linked to the Polyakov loopΦ = 1 N c (cid:42) Tr c P exp (cid:32) i (cid:90) β dτ A (cid:33)(cid:43) , (1)¯Φ = 1 N c (cid:42) Tr c P exp (cid:32) − i (cid:90) β dτ A (cid:33)(cid:43) . (2) Moreover, we include a coupling of the quarks to a mas-sive vector field ω . The resulting Lagrangian of the modelreads L =¯ q ( iγ µ D µ − g ( σ + iγ (cid:126)τ(cid:126)π ) − g ω γ µ ω µ ) q + 12 ( ∂ µ σ ) + 12 ( ∂ µ (cid:126)π ) − U m ( σ, (cid:126)π ) − U (Φ , ¯Φ; T ) − m ω ω + 14 F µν F µν , (3)where D µ = ∂ µ − iA µ , and F µν = ∂ µ ω ν − ∂ ν ω µ . Theparameters of the mesonic potential U m ( σ, (cid:126)π ) = λ σ + (cid:126)π ) + m σ + (cid:126)π ) − Hσ, (4)are tuned to vacuum properties of the σ and (cid:126)π mesons(see Appendix B). We use the Polyakov-loop potential U (Φ , ¯Φ; T ) determined in [61], by fitting quenched latticeQCD results for the equation of state and the Polyakov-loop susceptibilities. The parametrization of U (Φ , ¯Φ; T )is given in Appendix B.We compute the thermodynamic properties of thismodel including fluctuations of the scalar and pseu-doscalar meson fields within the framework of the FRGmethod. The ω meson, on the other hand, is treatedon the mean-field level. The MF treatment of the vec-tor field is justified by recent FRG results obtained atvanishing chemical potential indicating that the vectormeson mass remains above the cutoff during the FRGflow and that the temperature dependence of the screen-ing mass is very weak [62, 63]. Therefore, it is expectedthat fluctuations of the vector fields decouple from theflow. The Polyakov loop is also treated on the mean-fieldlevel. The Polyakov-loop variables and the vector fieldare tuned such that at the end of the calculation a sta-tionary point of the thermodynamic potential is reached.In the FRG framework, the effective average actionΓ k , which interpolates between the classical and the fullquantum action, is obtained by solving the renomaliza-tion group flow equation [56] ∂ k Γ k [ φ ] = 12 STr (cid:20)(cid:16) Γ (2) k [ φ ] + R k (cid:17) − ∂ k R k (cid:21) , (5)where φ denotes the quantum fields considered, STr is atrace over the fields, over momentum and over all internalindices. It also adds fermion contribution to the bosoncontribution with correct sign and prefactor. The reg-ulator function R k suppresses fluctuations at momentabelow k . Thus, effects of fluctuations of quantum fieldsare included shell by shell in momentum space, start-ing from a UV cutoff scale Λ. We employ the optimizedregulator introduced by Litim [64, 65], which yields analgebraic expression for the right hand side of the flowequation, Eq. (5). A. Thermodynamics at vanishing vector coupling
For vanishing vector coupling, the vector fields obvi-ously decouple completely. Hence, in such models thethermodynamics is controlled by the quark, as well as,the scalar and pseudoscalar fields [28, 29, 33, 66, 67].Using Ω k = Γ k /V | H =0 we obtain the thermodynamic po-tential densityΩ( T, µ q ) = − P ( T, µ q ) = Ω k =0 ( T, µ q , σ ∗ , Φ ∗ , ¯Φ ∗ ) − Hσ ∗ , (6)where µ q = µ B / k =0 is obtained by solving the flow equation, ∂ k Ω k ( T, µ q ; σ, Φ , ¯Φ)= k π (cid:26) E π coth (cid:18) E π T (cid:19) + 1 E σ coth (cid:18) E σ T (cid:19) − N c N f E q (cid:0) − N ( T, µ q ; σ, Φ , ¯Φ) − N ( T, − µ q ; σ, ¯Φ , Φ) (cid:1)(cid:27) , (7)with E q = (cid:112) k + g s σ , E σ = (cid:112) k + Ω (cid:48)(cid:48) k and E π = (cid:112) k + Ω (cid:48) k /σ . Here the prime denotes a derivative withrespect to σ field, g s is the Yukawa coupling, while N c = 3is the number of colors and N f = 2 the number of flavors.The Polyakov-loop modified quark occupation numbersare given by N ( T, µ q ; σ, Φ , ¯Φ)= 1 + 2 ¯Φ e ( E q − µ q ) /T + Φ e E q − µ q ) /T e ( E q − µ q ) /T + 3Φ e E q − µ q ) /T + e E q − µ q ) /T . (8)The flow equation (7) is solved numerically startingfrom the initial conditions specified at the momentumscale k = ΛΩ k =Λ = U m ( σ, | H =0 + U (Φ , ¯Φ; T ) . (9)The thermodynamic potential Ω( T, µ ), shown in Eq. (6),is obtained by evaluating the fully evolved Ω k =0 at thestationary point, SP = ( σ ∗ , Φ ∗ , ¯Φ ∗ ), determined by ∂ Ω k =0 ∂σ (cid:12)(cid:12)(cid:12)(cid:12) SP = H, ∂ Ω k =0 ∂ Φ (cid:12)(cid:12)(cid:12)(cid:12) SP = 0 , ∂ Ω k =0 ∂ ¯Φ (cid:12)(cid:12)(cid:12)(cid:12) SP = 0 . (10)Any thermodynamic quantity can be obtained by tak-ing the appropriate derivatives of the thermodynamic po-tential (6). However, quantities that require higher orderderivatives are difficult to extract by numerical differen-tiation, owing to limited numerical precision. We findthat improved numerical stability is obtained when thederivatives up to second order, ∂ Ω k ∂ Φ , ∂ Ω k ∂ ¯Φ , ∂ Ω k ∂ Φ , ∂ Ω k ∂ ¯Φ , ∂ Ω k ∂ Φ ∂ ¯Φ ,∂ Ω k ∂µ∂ Φ , ∂ Ω k ∂µ∂ ¯Φ , ∂ Ω k ∂µ , ∂ Ω k ∂µ , (11) are computed directly using flow equations, summarizedin Appendix A.We solve the coupled flow equations numerically bydiscretizing them on a grid of the order parameter σ . Inorder to recover the Stefan-Boltzmann limit, it is nec-essary to include the effect of fermion fluctuations abovethe UV cutoff Λ. This is done by adding the contributionof the thermal Polyakov-loop suppressed massless quarkloop, integrated from the UV cutoff to infinity [29]. Theinitial conditions for the flow equations are fixed by vac-uum criteria, as discussed in Appendix B.It should be noted that we also performed the RG cal-culations by expanding the potential in terms of Cheby-shev polynomials [68], instead of working on a uniformlydiscretized grid. This method is illustrated in Refs. [69–71]. Using the same initial conditions it yields the sameresults for all thermodynamic quantities, but up to ahigher numerical precision. This method was requiredto extract cumulants of high order. B. Thermodynamics at nonvanishing vectorcoupling
For nonvanishing vector coupling, the ω field is coupledto the chiral sector through the fermionic part of theflow. The temporal component of the ω field gains anonvanishing expectation value at nonzero density, whilethe spatial components vanish in a system with zero netquark current.An inspection of the Lagrangian (3), reveals that anonvanishing ω field effectively acts as a shift of the quarkchemical potential. Thus, it is convenient to introducean effective quark chemical potential, ν = µ − g ω ω [72],and continue using the formalism for vanishing vectorcoupling, albeit with modified initial conditions for theflow equationΩ k =Λ = U m ( σ, | H =0 + U (Φ , ¯Φ; T ) − m ω ω . (12)Since, for a given value of ω , this amounts to constantshift of the thermodynamic potential, it will not modifythe RG flow. Thus, the only effect of a nonzero ω field onthe flow is the replacement of the chemical potential µ by the previously defined effective chemical potential ν .The expectation value of ω is obtained by extremizingΩ k =0 . Since the flow for vanishing and nonvanishing vec-tor coupling can be related to each other, one can expressall results in terms of quantities computed at vanishingvector coupling. We denote the thermodynamic poten-tial density and net baryon density at vanishing vectorcoupling by ˜Ω( T, ν ) and ˜ n ( T, ν ), respectively. The ex-pectation value of the ω field is then determined by ∂ Ω k =0 ∂ω (cid:12)(cid:12)(cid:12)(cid:12) SP = 0 , = ⇒ g ω ω = − G ω ∂ ˜Ω ∂ν = G ω ˜ n ( T, ν ) , (13)where G ω = g ω /m ω . This yields the relation between thereal and the effective chemical potentials µ = ν + G ω ˜ n ( T, ν ) . (14)Moreover, the corresponding thermodynamic potentialfor nonzero vector coupling is given byΩ( T, µ ) = − P ( T, µ ) = ˜Ω(
T, ν ) − ( µ − ν ) G ω (15)We note, that for a first order phase transition, thisprocedure cannot be used to identify the critical value ofthermal parameters, owing to a flattening of the poten-tial, as a function of the order parameter. C. Net-baryon-number cumulants
Fluctuations of conserved charges, in particular of thenet baryon number, are valuable probes of critical phe-nomena and can be used to identify the location of thephase boundary of chiral symmetry restoration. Thefluctuations are characterized by cumulants, that can beobtained directly from the partition function. The n th -order baryon χ nB , and quark χ nq number cumulants are3 n χ nB = χ nq = 1 T − n ∂ n Ω( T, µ q ) ∂µ nq . (16)At vanishing vector coupling, the first order cumulant χ B , the net baryon density, can be obtained directly fromthe flow of ∂ Ω k /∂µ , while the second order cumulant χ B is given by χ B = 19 T (cid:18) ∂ Ω k =0 ∂µ q − (cid:88) i,j ∂ Ω k =0 ∂µ q ∂ϕ i (cid:18) ∂ Ω k =0 ∂ϕ∂ϕ (cid:19) − ij ∂ Ω k =0 ∂µ q ∂ϕ j (cid:19) , (17)where ϕ = ( σ, Φ , ¯Φ). The derivatives on the right handside of Eq. (17) are determined by the solution of flowequations introduced in Appendix A. Higher order cumu-lants are then computed using numerical differentiationof χ B .For nonvanishing vector coupling, the baryon numbercumulants can be expressed in terms of the correspondingcumulants ˜ χ nB for vanishing vector coupling. This is doneusing Eqs. (14) and (15). Based on Eq. (14) one can write dνdµ = 11 + 9 T G ω ˜ χ B ( T, ν ) , (18) which using Eq. (15), leads to the following relations χ B = ˜ χ B ,χ B = ˜ χ B T G ω ˜ χ B ( T, ν ) ,χ B = ˜ χ B (1 + 9 T G ω ˜ χ B ) ,χ B = ˜ χ B (1 + 9 T G ω ˜ χ B ) − T G ω (cid:0) ˜ χ B (cid:1) (1 + 9 T G ω ˜ χ B ) . (19)Here the cumulants ˜ χ kB are functions of temperature andthe reduced chemical potential ν . Using these relations,we can compute the influence of the vector interactionon the baryon number cumulants, using the RG flow ob-tained for vanishing vector coupling, G ω = 0. III. NET-BARYON-NUMBER CUMULANTSAND THE PHASE BOUNDARY
In the previous section we have introduced the methodemployed to compute cumulants of net baryon number inan effective chiral model in the presence of vector inter-actions. In the following, we discuss the properties ofthese cumulants near, and at the chiral phase boundary,at finite temperature and chemical potential.The chiral Lagrangian introduced above, shares thechiral critical properties with QCD. In particular, atmoderate values of the chemical potential, the PQMmodel exhibits a chiral transition belonging to the O (4)universality class [5, 12]. For larger values of µ , it revealsa Z (2) critical endpoint, followed by a first order phasetransition [12]. Consequently, the PQM model embodiesthe generic phase structure expected for QCD, with theuniversal O (4) and Z (2) criticality encoded in the scal-ing functions. Furthermore, due to the coupling of thequarks to the background gluon fields, the PQM modelincorporates ”statistical confinement”, i.e., the suppres-sion of quark and diquark degrees of freedom in thelow temperature, chirally broken phase [18, 30]. Con-sequently, by studying fluctuations of conserved chargesin the PQM model, one can explore the influence of chi-ral symmetry restoration and of ”statistical confinement”on the cumulants in different sections of the chiral phaseboundary. This study is of particular interest in the con-text of heavy-ion collisions, where cumulants of conservedcharges are expected to provide a characteristic signaturefor the QCD phase boundary and for the conjectured crit-ical endpoint [42–46, 49, 51]. Χ B6 (cid:144) Χ B2 g Ω = Χ B6 (cid:144) Χ B2 g Ω = g s (cid:144) Χ B4 (cid:144) Χ B2 g Ω = Χ B4 (cid:144) Χ B2 g Ω = g s (cid:144) s2 f Π Χ Ch - T (cid:144) T c FIG. 1. The temperature dependence of ratios of net-baryon-number cumulants, χ B /χ B and χ B /χ B in the PQM modelfor vanishing and nonzero vector coupling g ω . Also shown isthe chiral susceptibility, χ ch ( T ). The vertical line indicatesthe location of the chiral crossover temperature. A. Criticality of net-baryon-number cummulants inthe ( T, µ ) plane Generalized susceptibilities of conserved charges,which are given by higher-order derivatives of the ther-modynamic pressure with respect to the correspondingchemical potentials, may exhibit a nonmonotonic depen-dence on the thermodynamic parameters. This is par-ticularly the case in the vicinity of phase boundary andthe CEP. In the critical region of the chiral transition,the strength of the fluctuations and the sign of the sus-ceptibilities are by and large determined by the singularpart of the free energy, which is encoded in the universalscaling functions, common to QCD and the PQM model.Thus, generic structures of the susceptibilities and rela-tions between them near the phase boundary, can alsobe studied in the PQM model. Of particular interest isthe behavior of susceptibilities along the O (4) crossoverline, and their modification as the CEP is approached.In the model calculations, the position of the CEP de-pends on the strength of the vector interaction. Thus,by changing the vector coupling g ω , one can assess thedependence of critical fluctuations in different sectors ofthe phase boundary on the location of the CEP.For nonvanishing light-quark masses, the chiral sym-metry is explicitly broken, which implies that at finitetemperatures and moderate values of the baryon chemi-cal potential, the system exhibits a chiral crossover tran-sition. Thus, at small µ q , the fluctuations of conservedcharges remain finite. However, because QCD matterat the physical values of the u and d quark masses iswithin the critical region of the second order phase transi-tion, the fluctuations are still influenced by O (4) critical-ity [6, 73–76]. At physical quark masses, a genuine phasetransition in QCD and the associated singular behavior of fluctuations are only expected at the conjectured CEPand along the line of the first-order phase transition.To explore the phase diagram with net-baryon-numberfluctuations, we need to know the qualitative dependenceof net-baryon-number cumulants on the temperature andbaryon chemical potential. In Fig. 1 we show the T -dependence at µ q = 0, and in Figs. 2 and 3, the contourplots in the ( T, µ q )-plane of ratios of net-baryon-numbersusceptibilities. We focus on the following ratios of net-baryon-number susceptibilities: χ n,mB = χ nB ( T, µ B ) χ mB ( T, µ B ) , (20) κσ = χ B ( T, µ B ) χ B ( T, µ B ) , S B σ M = χ B ( T, µ B ) χ B ( T, µ B ) , (21)where M is the mean, σ the variance, S B the skewnessand κ the kurtosis of the net-baryon-number distribution.The statistical descriptors are related to the susceptibil-ities through M = V T χ B , σ = V T χ B etc..At vanishing chemical potential, all odd susceptibili-ties of net baryon number vanish, owing to the baryon-antibaryon symmetry. In addition, in the O (4) universal-ity class, the second and fourth order cumulants remainfinite at the phase transition temperature at µ q = 0 evenin the chiral limit, implying that only sixth and higherorder susceptibilities diverge. Thus, for physical quarkmasses, only higher order cumulants, χ nB with n > O (4) criticality at µ q = 0 [46]. A furtherconsequence of the baryon-antibaryon symmetry is theequality of the ratios χ m − , n − B = χ m, nB (22)for any integer m and n ≥ µ q = 0. For χ , B χ , B ,the equality at small µ q can also be confirmed by a directcomparison of the right panel of Fig. 2 with the left panelof Fig. 3.At finite net baryon density the singularity at the O (4)critical line is stronger than at µ q = 0. Thus, in thiscase the third-order cumulant and all higher-order ones,diverge at the O (4) line. The second order cumulant χ B remains finite, and only diverges at the tricritical pointfor vanishing quark masses, and at the CEP for nonzeroquark masses.In Fig. 1 we show results on the temperature depen-dence of several ratios χ n,mB of net-baryon-number sus-ceptibilities, together with the variance of the chiral con-densate, χ ch . The location of the maximum of the chiralsusceptibility, χ ch , defines the pseudo-critical tempera-ture, T c .As shown in Fig. 1, the kurtosis κσ exhibits a rapiddrop by approximately an order of magnitude at thephase boundary. This strong reduction of the kurtosisis attributed to the deconfinement transition, where thedegrees of freedom carrying baryon number change frombaryons to quarks [49, 50]. A qualitative understandningof this issue is obtained in the Boltzmann approximation, Phase Boundary Μ q (cid:144) T c T (cid:144) T c Χ B1 (cid:144) Χ B2 , g Ω = - Phase Boundary Μ q (cid:144) T c T (cid:144) T c Χ B3 (cid:144) Χ B1 , g Ω = FIG. 2. Contour plots of the ratios χ B /χ B and χ B /χ B in the ( T, µ )-plane, computed in the PQM model. The broken linesindicate the location of the chiral crossover phase boundary. - Phase Boundary Μ q (cid:144) T c T (cid:144) T c Χ B4 (cid:144) Χ B2 , g Ω = - Phase Boundary Μ q (cid:144) T c T (cid:144) T c Χ B4 (cid:144) Χ B2 , g Ω = s FIG. 3. Contour plots of the kurtosis ratio, χ B /χ B , in the ( T, µ q )-plane. The left-hand figure corresponds to PQM modelresults obtained with vanishing vector coupling g ω = 0, while the right-hand one shows results obtained with g ω = 0 . g s . where the baryon contribution to the thermodynamicpressure is of the form P B ( T, µ B ) (cid:39) F ( T ) cosh( Bµ B /T ).Here B is the baryon number of the relevant degrees offreedom and F ( T ) a function of the corresponding ex-citation spectrum. At low temperatures, the degrees offreedom are baryons, with B = ±
1, while at high tem-peratures they are quarks with B = ± /
3. Consequently,the kurtosis ratio, κσ = χ B /χ B , is approximately pro-portional to the square of the baryon number of the rele-vant degrees of freedom [49, 50]. This clarifies the causefor the rapid change of κσ at the crossover transitionfrom (cid:39) (cid:39) /
9, seen in Fig. 1. In the limit where thecharacteristic mass is small compared to the tempera-ture, the kurtosis ratio is modified by quantum statistics,to κσ = (1 / /π ). We note, that the kurtosis ratiois independent of the mass spectrum as well as of anykinematic cuts, as long as the Boltzmann limit remainsvalid. The ratio χ , B , also shown in Fig. 1, exhibits morestructure at the chiral transition. The observed char-acteristic temperature dependence of this ratio, in par-ticular with a region of negative values at T (cid:39) T c , is aconsequence of the residual chiral O (4) criticality [46].In the absence of chiral critical fluctuations, χ , B wouldshow a similar behavior as κσ , with a smooth reductionfrom unity at low temperatures towards zero above thedeconfinemnt transition. The distinctive behavior of the χ , B and χ , B ratios was already obtained in the PQMmodel within the FRG approach [46], and agrees quali-tatively with LQCD results [77]. Thus, the PQM modelcorrectly captures the physics of QCD related to decon-finement and to the critical dynamics at the chiral tran-sition. The negative values of χ , B at the chiral crossoverwere proposed as a signature for partial restoration ofchiral symmetry in heavy ion collisions [46].As shown in Fig. 1, in the presence of repulsive inter-actions, χ , B and the kurtosis ratios are modified. Thevector interaction leads to a downward shift in tempera-ture of the ratios relative to the phase boundary, as wellas, a suppression of the sixth order susceptibility in thetemperature range below T c . Nevertheless, the qualita-tive form of the ratios is preserved. In particular, thecharacteristic structure, where the sixth order cumulantis negative in a range of temperatures near T c owing to O (4) criticality, is not eliminated by the repulsive inter-action.In Fig. 2 we show contour plots of the ratios χ , B and χ , B in the ( T, µ q )-plane. As noted above, all odd cumu-lants vanish at µ q = 0, owing to baryon-antibaryon sym-metry. Consequently, χ , B | µ q =0 = 0 for any T , while theratio χ , B | µ q =0 is nonvanishing. At low and high T , rela-tive to T c , this ratio is consistent with unity and 2 / (3 π ),respectively, as expected since χ , B | µ q =0 = χ , B | µ q =0 . Asindicated in Figs. 1 and 2, the ratio χ , B decreases withtemperature, and depends weakly on the chemical po-tential. Moreover, Fig. 2, shows that for µ q < T , χ , B increases with µ q , and is only weakly dependent on thetemperature. Thus, the ratio χ , B can be used as a mea-sure of the temperature, while χ , B provides a gauge ofthe chemical potential.At small µ q /T , the properties of the first four suscep-tibilities, χ nB with n = 1 , ..,
4, and consequently theirratios near the chiral crossover are dominantly affectedby the coupling of the quarks to the Polyakov loop, andthe resulting statistical confinement. The critical chiraldynamics, i.e. the O (4) and Z (2) criticality at the chiralcrossover transition and at the CEP, respectively, unfoldsat larger µ q /T . Near the CEP, there is a strong varia-tion of the cumulants with T and µ q , which increaseswith the order of cumulants. As shown in Figs. 2 and 3,the qualitative behavior of the cumulant ratios on linesgoing to the CEP depends strongly on the direction fromwhich the CEP is approached. This behavior is governedby the critical properties encoded in the O (4) and Z (2)universal scaling functions.The influence of the CEP on the characteristics of thevarious cumulant ratios can be studied by varying thevalue of the vector coupling, thus changing the positionof the CEP in the ( T, µ q )-plane. This is illustrated inFig. 3 for the kurtosis ratio. A comparison of the resultsshown in the left and right contour plots shows that withincreasing g ω , the curvature of the phase boundary is re-duced and the position of the CEP is shifted to lower T and larger µ q [78]. As is clearly seen, when comparing theleft and right plots of Fig. 3, there is a corresponding shiftof the contours of χ , B . These results clearly illustratethe influence of the critical endpoint on the fluctuationobservables. A shift of the CEP to larger net baryon den-sity, suppresses the magnitude of the net-baryon-numbersusceptibilities at a given T and µ q . B. Net-baryon-cumulant ratios and freeze-out inheavy ion collisions
In heavy-ion collisions, the thermal fireball formed inthe quark-gluon plasma phase undergoes expansion andpasses through the QCD phase boundary at some point( T c , µ c ) point, which depends on the collision energy, √ s . Analysis of ratios of particle multiplicities indi-cate that at high beam energies (small values of µ q /T ),the freeze-out occurs at, or just below the phase bound-ary. Thus, the beam energy dependence of net-baryon-number susceptibilities can provide insight into the struc-ture of the QCD phase diagram and information on theexistence and location of the CEP. Consequently, it is ofphenomenological interest to compute the properties offluctuations of conserved charges along the chiral phaseboundary. Since there, the critical structure and the re-lations between different susceptibilities are governed bythe universal scaling functions, the generic behavior ofratios of net-baryon-number susceptibilities can be ex-plored also in model calculations.In Fig. 4 we show χ , B , as well as skewness and kurtosisratios of net baryon number, obtained in the PQM modelalong the chiral phase boundary, with and without vectorrepulsion. In the former case, the position of the CEP isshifted to smaller T and larger µ q .The ratio χ , B exhibits a maximum along the phaseboundary. For small ( µ q /T ) c this ratio is well approxi-mated by tanh(3 µ q /T ), while after reaching a maximum,it decreases as the CEP is approached. The cumulants χ B and χ B remain finite along the O (4) line. However, atthe CEP the variance of the net-baryon-number fluctua-tions diverges. This is the cause for the observed decreaseof χ , B at larger ( µ q /T ) c . For nonzero vector repulsion,this reduction is weakened, owing to the shift of the CEPto lower temperature. We note, that for ( µ q /T ) c < . χ , B ratio is hardly modified by the repulsive interac-tions. This can be traced back to a cancellation betweenthe suppression of χ B at nonzero g ω (19) and the shiftof the chemical potential (14); both axes are scaled bythe factor (1 + 9 T G ω ˜ χ B ( T, ν )), for small µ q . On theother hand, the ratio χ B /χ B is reduced by the factor(1 + 9 T G ω ˜ χ B ( T, ν )), for nonvanishing vector interac-tion. Hence, this ratio is reduced at all values of µ q /T .As seen in Fig. 4, the skewness ratio, χ , B , is an in-creasing function of ( µ q /T ) c along the phase boundaryand at the CEP it diverges [79]. There is also a clearsuppression of this ratio along the phase boundary dueto repulsive interactions, as seen also for χ , B . At small( µ q /T ) c = 0, the skewness and kurtosis ratios are equalto each other, and differ very little up to ( µ q /T ) c (cid:39) . µ q /T ) c , the kurtosis is decreasing and skew-ness increasing along the phase boundary, in agreementwith recent LGT results [80]. This behavior also reflectstheir properties near the critical point. As the CEP isapproached along the phase boundary, the kurtosis di-verges to minus infinity, while the skewness diverges to Χ B3 (cid:144) Χ B2 g Ω = Χ B3 (cid:144) Χ B2 g Ω = g s (cid:144) Χ B1 (cid:144) Χ B2 g Ω = Χ B1 (cid:144) Χ B2 g Ω = g s (cid:144) Μ q (cid:144) T Χ B4 (cid:144) Χ B2 g Ω = Χ B4 (cid:144) Χ B2 g Ω = g s (cid:144) Χ B3 (cid:144) Χ B1 g Ω = Χ B3 (cid:144) Χ B1 g Ω = g s (cid:144) - - Μ q (cid:144) T FIG. 4. Ratios of cumulants of net-baryon-number fluctuations in the PQM model computed on the chiral phase boundary withand without vector repuslion. The solid and long-dashed lines correspond to g ω = 0, and the dash-dotted lines to g ω = g s / plus infinity, as noted above. Thus, both ratios becomesless singular as g ω is increased, i.e., as the CEP is shiftedto larger µ q .The characteristics of the various net-baryon-numbersusceptibilities on the phase boundary, shown in Fig. 4,are expected to be similar in QCD. This is because, theyare, by and large, determined by O (4) critical fluctua-tions and by ”confinement”, which are both common toQCD and the PQM model. This opens the possibilityto verify directly, if these features of criticality are alsoreflected in the data on net-proton-number fluctuationsobtained in heavy-ion collisions by the STAR Collabora-tion at several RHIC energies [52–54].Clearly, a direct comparison of model results with datahas to be taken with caution. Although, the model pro-vides a viable description of the dynamics that drives thesystem towards chiral symmetry restoration, the spec-trum of hadronic degrees of freedom in the low tempera-ture phase is incomplete. Moreover, net-baryon-numberfluctuations are in nucleus-nucleus collision experimentsquantified by the net proton number. It has been exten-sively discussed, to what extent are net-proton-numberfluctuations accurate proxyies for those of the net baryonnumber [81–84]. Furthermore, there are kinematical cuts,imposed on the STAR data on the cumulants of net pro-ton number, which are not accounted for in the modelresults. However, these differences are, to a large extent,eliminated by considering ratios of susceptibilities. Thisassumption is supported by the behavior of the ratio χ , B ,which, in spite of the differences discussed above, is wellapproximated by tanh(3 µ q /T ) (see Fig. 5) in LQCD, inthe PQM model as well as in the STAR data. The factthat χ , B is well approximated by this functional formindicates that in the transition region, the effective de-grees of freedom with nonvanishing baryon number have B = ± √ s and the thermal parameters ( µ q , T ).Here we employ the phenomenological relation, obtainedby analysing the freeze-out conditions in terms of thehadron-resonance-gas model (HRG) [38–41]. We then usethe resulting dependence of µ B and T on √ s to assign avalue for the ratio ( µ q /T ) to each of the STAR beam en-ergies. We note that, for µ q /T <
1, the phenomenologicalfreeze-out line coincides within errors with the crossoverphase boundary obtained in lattice QCD [3, 40]. Thismotivates a comparison of model results on net-baryon-number fluctuations near the phase boundary with data.Such an analysis was first done using LQCD results inRef. [60].In Fig. 5, we show the STAR data on net-proton-number susceptibility ratios and the corresponding PQMmodel results on net-baryon-number fluctuations com-puted along the phase boundary. The model results forthe ratios χ , B , χ , B and χ , B are in qualitative agreementwith the data in the whole energy range. For the kurtosisratio, χ , B , this is the case also up to the SPS energy, i.e.,for √ s ≥
20 GeV. However, for µ q /T > .
5, the data onthe kurtosis ratio exhibits a qualitatively different depen-dence on µ q /T than expected for the critical behavior of χ , B , as the CEP is approached along the phase bound-ary.As noted above, the ratio χ , B is, on the phase bound-ary, approximately given by χ , B (cid:39) tanh(3 µ q /T ) up to µ q /T (cid:39)
1, as seen in Fig. 5. This form of the ratio ofthe lowest order cumulants is also obtained in the HRGmodel and is consistent with LQCD [48, 51]. Fig. 5 re-veals that the data are consistent with this form as well.This fact clearly supports the use of the √ s dependenceof µ q /T obtained from the HRG model analysis of mul-tiplicities in the comparison of the PQM model resultswith data.In order to assess the sensitivity of the fluctuation ob-servables to model parameters, we show in Fig. 5 ratios Χ B1 (cid:144) Χ B2 Χ B3 (cid:144) Χ B2 Set ASet B Set CSet DTanh @ Μ B (cid:144) T D Μ q (cid:144) T s @ GeV D Χ B4 (cid:144) Χ B2 Set ASet B Set CSet D Χ B3 (cid:144) Χ B1 - Μ q (cid:144) T s @ GeV D FIG. 5. Ratios of cumulants of net-baryon-number fluctuations in the PQM model, computed along the chiral phase boundary,for four sets of model parameters (for details see Appendix B). Also shown are the preliminary STAR data [53, 54], assumingthe relation between the ratio ( µ q /T ) and the collision energy obtained by analysing the chemical freeze-out conditions [39–41].The green full line in the left-hand figure shows the baseline result, χ , B = tanh(3 µ q /T ). Χ B1 (cid:144) Χ B2 Χ B3 (cid:144) Χ B2 Set ASet B Set CSet D Μ q (cid:144) T s @ GeV D Χ B4 (cid:144) Χ B2 Set ASet B Set CSet D Χ B3 (cid:144) Χ B1 - Μ q (cid:144) T s @ GeV D FIG. 6. Ratios of cumulants of net-baryon-number fluctuations in the PQM model along the freeze-out line, obtained byfitting χ B /χ B to the STAR data. The four sets of model parameters used and the preliminary STAR data shown, are the sameas in Fig. 5. of cumulants for four sets of model parameters. The dif-ferent sets are obtained by varying the sigma meson massand the form of the Polyakov-loop potential. The param-eter sets are described in Appendix B. Fig. 5 shows that,although some quantitative differences can be identifiedat larger µ q /T , the skewness, kurtosis and χ , B ratiosalong the phase boundary are qualitatively similar forthe different sets of model parameters.In the comparison of model predictions with data inFig. 5, we assume, that the freeze-out of the net-baryon-number fluctuations, tracks the chiral phase boundary.Clearly, this simple assumption provides a qualitative un- derstanding of the data. In order to obtain a more quan-titative description, we follow Refs. [60] and [85, 86], anddetermine the freeze-out conditions by fitting the dataon the χ , B ratio, using the ( √ s )-dependence of µ q /T obtained from the fit of the HRG model to particle mul-tiplicities [39–41].In Fig. 6 we show the fluctuation ratios along thefreeze-out line, which is fixed through the skewness data.The model results are obtained for the four sets of initialconditions. Fig. 6 clearly shows that, along the freeze-out line, the spread of all fluctuations ratios consideredfor the various parameter sets is much weaker than that0 Χ B6 (cid:144) Χ B2 Χ B5 (cid:144) Χ B1 Χ B4 (cid:144) Χ B2 s2 f Π Χ Ch - - T (cid:144) T Χ B6 (cid:144) Χ B2 Χ B5 (cid:144) Χ B1 Χ B4 (cid:144) Χ B2 Χ B3 (cid:144) Χ B1 - - Μ q (cid:144) T s @ GeV D FIG. 7. Left-hand figure: The temperature dependence of ratios of net-baryon-number cumulants and the chiral susceptibility χ ch , all computed in the PQM model at µ q = 0 . χ B /χ B , χ B /χ B and χ B /χ B along the freeze-out line obtained by fitting χ B /χ B (also shown)to data. These results were computed using the model parameters of set A (see Appendix B). The bands about χ B /χ B and χ B /χ B reflect the experimental error in χ B /χ B , which leads to an uncertainty in the freeze-out temperature. observed in Fig. 5 along the phase boundary. This indi-cates, that moderate changes of the sigma mass and mod-ifications of the form of the Polyakov loop potential maylead to a shift in the temperature scale but essentiallywith no change of the relative structure of the cumulantratios.The results presented in Fig. 6 clearly show thatthe model provides a very good description of the dataon χ , B and χ , B . Also the kurtosis data, obtained athigher collision energies, are consistent with model re-sults. However, at √ s <
20 GeV they exhibit a differenttrend, with the data increasing rapidly at lower energies,while the model result keeps decreasing. We concludethat an increase of χ , B ratio beyond unity, observed inthe STAR data at √ s <
20 GeV, is not expected in equi-librium on the chiral critical line nor on the freeze-outline.As noted above, the ratio χ , B is particularly interest-ing for identifying criticality governed by the O (4) univer-sality class. This is seen in Fig. 1 at µ q = 0, where χ , B isnegative at T = T c . At µ q >
0, the influence of criticalityis even more pronounced, as shown in Fig. 7. There, χ , B exhibits a highly nonmonotonic structure near the phaseboundary.In the left panel of Fig. 7 we show the temperaturedependence of the ratio χ , B near the phase boundary for µ q = 0 . T c and develops a maximum just above T c . Thevalue of χ , B is thus very sensitive to the freeze-out tem-perature.On the right of Fig. 7 we show the ratios χ B /χ B , χ B /χ B and χ B /χ B computed, using the model param- eters of set A, along the freeze-out line which is de-termined by fitting χ B /χ B to data. At vanishing µ q , χ B /χ B = χ B /χ B and χ B /χ B = χ B /χ B , as a conse-quence of Eq. (22), while at larger µ q /T , these ratiosseparate.We note at this point, that the equality of the ratios χ B /χ B and χ B /χ B in the STAR data at the highest en-ergy is a strong indication that the fluctuations probedby these cumulants are in thermal equilibrium. It is veryunlikely that a system not in equilibrium would yieldratios of cumulants that satisfy Eq. (22). Note, that at µ q = 0, the critical O (4) fluctuations yield divergent con-tributions only to χ B and higher cumulants [46]. Thus,the fluctuations probed by χ B /χ B and χ B /χ B are notcritical. However, a measurement of the ratios involvingthe fifth and sixth order cumulants would probe whetherthe O (4) critical fluctuations are in equilibrium or not.Obviously, this test of equilibration is meaningful only atsmall µ q /T , i.e., only at the highest energies.At moderate values of µ q /T , the χ B /χ B ratio is neg-ative and deviates clearly from χ B /χ B . At still lowerenergies, it exhibits a strong increase towards the CEP,where it diverges. Similarly, χ B /χ B decreases at mod-erate µ q /T , and increases strongly as the CEP is ap-proached. These results indicate that in heavy ion col-lisions χ B /χ B and χ B /χ B will exhibit strong nonmono-tonic dependencies on √ s .Recently, first results on χ B /χ B in Au-Au collisions at √ s = 200 GeV where reported by the STAR Collabora-tion for several centralities [87]. The data show a strongsuppression of fluctuations compared to the kurtosis ra-tio. In mid-central and the most central collisions, the χ B /χ B fluctuation ratio is negative, albeit with still verylarge statistical uncertainties. In fact, given the large er-1rors for the most central collisions, the preliminary datais consistent with a vanishing χ B . A value close to zerois consistent with the model results, shown in the rightpanel of Fig. 7.The comparison of model results on ratios of net-baryon-number susceptibilities with the STAR data inFigs. 5, 6, and 7, shows that the data, with the excep-tion of kurtosis at low energies, follow general trends ex-pected due to critical chiral dynamics. We note, thatthe ratios of net-baryon-number susceptibilities near thephase boundary involving net-baryon number cumulants χ nB with n ≥ O (4) and Z (2) universality classes, respec-tively. This observation indicates, that by measuringfluctuations of conserved charges in heavy-ion collisions,we are indeed probing the QCD phase boundary, andthus accumulating evidence for chiral symmetry restora-tion.However, as discussed above, there are several uncer-tainties and assumptions which must be thoroughly un-derstood before the QCD phase boundary can be pinneddown with confidence. Possible contributions to fluctu-ation observables from effects not related with criticalphenomena, like e.g. baryon-number conservation [82]and volume fluctuations [88–93] are being explored. Wenote, however, that volume fluctuations do not modifythe double ratios proposed as tests of equilibration at µ q →
0. We also mention the rather strong sensitivityof higher order net-proton-number cummulants on thetransverse momentum range imposed in the analysis ofthe STAR data. Nevertheless, it is intriguing, that thedynamics of the model, provides a good description ofthe STAR data (except for χ B at the lowest energies),without all these effects of non-critical origin. It remainsan important task to assess the effect of theses additionalsources of fluctuations in the whole energy range probedby the expriments. IV. SUMMARY AND CONCLUSSIONS
We have studied the influence of the chiral phase tran-sition on fluctuation observables in strongly interactingmedium at finite temperature and net baryon numberdensity. We focused on the properties of net-baryon-number fluctuations which are quantified by the n th -order susceptibilities, χ nB . The cumulants χ nB are directlyinfluenced by the chiral phase transition due to the cou-pling of the quarks to the scalar sigma field. Further-more, since the χ nB are accessible experimentally, theyare ideal observables to identify the phase boundary andthe critical structures in the QCD phase diagram.The dynamics was modelled with the Polyakov loopextended Quark-Meson (PQM) model. To correctly ac-count for the critical behavior at the chiral symmetryrestoration transition in the O (4) and Z (2) universal-ity classes, we employed the Functional RenormalizationGroup (FRG). We formulated the FRG equations in the presence of repulsive interactions, and derived the flowequations for derivatives of the thermodynamic pressure.The main point of our studies was to identify the re-lations between χ nB susceptibilities in different ( T, µ B ) -regions of the phase diagram. To reduce the influence ofthe non-critical characteristics of the model, like e.g. themass spectrum or the kinematical cuts on particle mo-mentum distributions, we have computed ratios of sus-ceptibilities. Here, of particular interest, are ratios ofthe first and second order cumulants, χ B /χ B , the skew-ness, χ B /χ B , the kurtosis, χ B /χ B , and the sixth-ordercumulant χ B /χ B . The higher-order ( n ≥
3) cumulantsare probes of chiral criticality and can be reconstructedexperimentally in nucleus-nucleus collisions.We have calculated, for the first time, the contour plotsin the (
T, µ B )-plane of the above ratios in the PQMmodel within FRG approach. We have quantified thesystematic relations along the phase boundary and thephenomenological freeze-out line, extracted by fitting theskewness data. The influence of the repulsive interactionsand the position of the critical endpoint on different χ nB ratios, was also explored for the first time.The model results were confronted with ratios of cumu-lants of net proton number measured by the STAR Col-laboration in nucleus-nucleus collisions at energies rang-ing from √ s = 7 . √ s -dependence of µ q /T , follows thephenomenological freeze-out conditions in nucleu-nucleuscollisions. The validity of this assumption is supportedby the observation that, for µ q /T <
1, the ratio of thefirst and second order cumulants is well approximatedby χ B /χ B = tanh(3 µ q /T ), in the model, in LQCD andin the data, using the adopted relation between beamenergy and µ q /T .We have shown, that the STAR data for the ratioscomposed of χ nB with n = 1 , , χ B /χ B to data. The data on the kurtosis ra-tio are also consistent with model systematics at ener-gies √ s ≥
20 GeV. However, the strong enhancement of κσ >
1, observed in nucleus-nucleus collisions at lowerenergies, is not reproduced by the model along the pathin the ( µ q /T )- plane, where the ratios of lower order cu-mulants, χ nB with n = 1 , ,
3, are well described. We have2shown, that this conclusion is not affected by the initialconditions nor by the parametrization of the Polyakov-loop potential.Clearly, the comparison of model results and data isbiased by the various assumptions and uncertainties, dis-cussed in Sec. III. However, in spite of all uncertainties,the ratios of net-proton-number susceptibilities obtainedby STAR are, as shown in this study, largely consis-tent with the systematics expected near the chiral phaseboundary.It is interesting to note that a recent analysis of theSTAR data on skewness and kurtosis ratios, at √ s ≥ . µ q /T ,shows that the data exhibits all features expected in QCDnear the phase boundary [80]. This result support ourfindings in PQM model that near chiral phase boundary, ratios of net-baryon-number cumulants, capture proper-ties, which are common to QCD and the PQM model. ACKNOWLEDGMENTS
We acknowledge stimulating discussions with PeterBraun-Munzinger, Frithjof Karsch and Nu Xu. G.A. isalso grateful to Vladimir Skokov for discussions on vec-tor interactions. The work of B.F. and K.R. was partlysupported by the Extreme Matter Institute EMMI. K.R. also acknowledges partial supports of the Polish Na-tional Science Center (NCN) under Maestro grant DEC-2013/10/A/ST2/00106. G. A. acknowledges the supportof the Hessian LOEWE initiative through the HelmholtzInternational Center for FAIR (HIC for FAIR).
Appendix A: Flow equations for derivatives of the thermodynamic potential
In the following, we introduce the flow equations that are used to calculate χ B and χ B cumulants of net-baryon-number fluctuations. We start from the flow equation (7) for the grand canonical potential, ddk Ω k = 3 F b (cid:0) m π (cid:1) + F b (cid:0) m σ (cid:1) + F f (cid:0) µ, Φ , ¯Φ (cid:1) , (A1)where the boson and fermion contribution is introduced, as F b ( m ) = k π coth (cid:16) √ k + m T (cid:17) √ k + m , (A2) F f ( µ, Φ , ¯Φ) = − k N c N f π E q (cid:0) − N ( T, µ ; σ, Φ , ¯Φ) − N ( T, − µ ; σ, ¯Φ , Φ) (cid:1) (A3)with the pion and sigma masses obtained from m π = 1 σ ∂ Ω k ∂σ , m σ = ∂ Ω k ∂σ . (A4)The flow equations for derivatives, introduced in Eq. (11), are derived for each gridpoint from Eq. (A1), as ddk ∂ Ω k ∂λ i = ∂∂λ i d Ω k dk = 3 F (cid:48) b (cid:0) m π (cid:1) σ ∂ Ω k ∂λ i ∂σ + F (cid:48) b (cid:0) m σ (cid:1) ∂ Ω k ∂λ i ∂ σ + ∂∂λ i F f ( µ, Φ , ¯Φ) ,ddk ∂ Ω k ∂λ i ∂λ j = ∂ ∂λ i ∂λ j d Ω k dk = 3 F (cid:48)(cid:48) b (cid:0) m π (cid:1) σ ∂ Ω k ∂λ i ∂σ ∂ Ω k ∂λ j ∂σ + 3 F (cid:48) b (cid:0) m π (cid:1) σ ∂ Ω k ∂λ i ∂λ j ∂σ + F (cid:48)(cid:48) b (cid:0) m σ (cid:1) ∂ Ω k ∂λ i ∂ σ ∂ Ω k ∂λ j ∂ σ + F (cid:48) b (cid:0) m σ (cid:1) ∂ Ω k ∂λ i ∂λ j ∂ σ + ∂ ∂λ i ∂λ j F f ( µ, Φ , ¯Φ) , (A5)where λ i and λ j stand for one of the { µ, Φ , ¯Φ } parameters, and the prime denotes derivative with respect to m .3 a a a a a -44.14 151.4 -90.0677 2.77173 3.56403 b b b b -0.32665 -82.9823 3.0 5.85559 c c c c c -50.7961 114.038 -89.4596 3.08718 6.72812 d d d d d Appendix B: Polyakov-loop potential and initialconditions
To verify possible model dependence of our results,we have considered different parameterizations of thePolyakov-loop potential.We apply the polynomial Polyakov-loop potential [19],as U (Φ , ¯Φ; T ) T = − b ( T )2 (cid:0) Φ ¯Φ (cid:1) − b (cid:0) Φ + ¯Φ (cid:1) + b (cid:0) Φ ¯Φ (cid:1) , (B1)where b = 0 . b = 7 .
5, and b ( T ) = a + a (cid:18) T T (cid:19) + a (cid:18) T T (cid:19) + a (cid:18) T T (cid:19) , (B2)with a = 6 . a = − . a = 2 . a = − .
44, and T = 270 MeV.The parameter T corresponds to the value of the crit-ical temperature at the first order deconfinement phasetransition in a pure Yang-Mills theory. Effects of un-quenching can be taken into account by changing T . Inthe case of two quark flavors, we use T = 208 MeV fromRef. [28].The second potential applied in our calculations, ac-counts for the color group structure, and reproduces lat-tice results on the equation of state and fluctuations ofthe Polyakov loops, in a pure SU(3) gauge theory. The potential is parameterized [61], as U (Φ , ¯Φ; T ) T = − a ( T ) (cid:0) Φ ¯Φ (cid:1) + b ( T ) log M H (Φ , ¯Φ)+ 12 c ( T ) (cid:0) Φ + ¯Φ (cid:1) + d ( T ) (cid:0) Φ ¯Φ (cid:1) , (B3)where M H (Φ , ¯Φ) = (cid:16) −
6Φ ¯Φ + 4 (cid:0) Φ + ¯Φ (cid:1) − (cid:0) Φ ¯Φ (cid:1) (cid:17) . (B4)The temperature dependent coefficients are given by a ( T ) = a t + a t + a t + a t + a , b ( T ) = b t − b (cid:16) − e b /t b (cid:17) ,c ( T ) = c t + c t + c t + c t + c , d ( T ) = d t + d t + d t + d t + d , (B5)where t = T /T . The numerical values of the con-stants in Eq. (B5), are summarized in Tab. I. We set T = 270 MeV for this potential. Lower values of T would result in distinct position of the peaks of thePolyakov loop and chiral susceptibilities, which is in con-trast with LQCD results.The initial conditions for the flow equations are deter-mined by the vacuum conditions: (cid:104) σ (cid:105) T =0 ,µ =0 = f π = 93 MeV ,m π = 138 MeV , m q = 300 MeV . (B6)The Yukawa coupling g s and the external field are thendetermined from g s = m q /f π , H = m π f π . (B7)However, the values of λ and m , in the initial meson po-tentials, are still undefined. They can be fixed by spec-ifying the cutoff scale Λ, where the flow is started, andthe mass of the sigma meson, m σ . We apply, two distinctvalues of the initial cutoff scales, and the sigma masses,as Λ = 700 MeV , (B8) m σ ≈
410 MeV , λ = − , m = 2 . · MeV , Λ = 950 MeV , (B9) m σ ≈
500 MeV , λ = 1 . , m = 5 . · MeV . where the parameters in Eq. (B9), were introduced inRef. [94].From the two parametrization of the Polyakov-loop po-tentials, Eqs. (B1) and (B3), and the two chiral initialconditions defined in Eqs. (B8) and (B9), we define inTable II, the four different parameter sets that are usedin our calculations to study the influence of the modelassumptions on the final results. If the parameter set isnot explicitly specified, then results refer to set A.4 [1] Y. Aoki, G. Endrodi, Z. Fodor, S. D. Katz, and K. K.Szabo, Nature , 675 (2006).[2] Y. Aoki, S. Borsanyi, S. Durr, Z. Fodor, S. D. Katz,S. Krieg, and K. K. Szabo, JHEP , 088 (2009).[3] S. Borsanyi, Z. Fodor, C. Hoelbling, S. D. Katz, S. Krieg,C. Ratti, and K. K. Szabo (Wuppertal-Budapest), JHEP , 073 (2010).[4] A. Bazavov et al. , Phys. Rev. D , 054503 (2012).[5] R. D. Pisarski and F. Wilczek, Phys. Rev. D , 338(1984).[6] S. Ejiri, F. Karsch, E. Laermann, C. Miao, S. Mukherjee, et al. , Phys. Rev. D , 094505 (2009).[7] M. Asakawa and K. Yazaki, Nucl. Phys. A , 668(1989).[8] A. Barducci, R. Casalbuoni, S. De Curtis, R. Gatto, andG. Pettini, Phys. Lett. B , 463 (1989).[9] F. Wilczek, Int. J. Mod. Phys. A , 3911 (1992).[10] A. M. Halasz, A. D. Jackson, R. E. Shrock, M. A.Stephanov, and J. J. M. Verbaarschot, Phys. Rev. D , 096007 (1998).[11] H.-T. Ding, F. Karsch, and S. Mukherjee, Int. J. Mod.Phys. E , 1530007 (2015).[12] B.-J. Schaefer and J. Wambach, Phys. Rev. D , 085015(2007).[13] Z. Fodor and S. D. Katz, Phys. Lett. B , 87 (2002).[14] C. R. Allton, S. Ejiri, S. J. Hands, O. Kaczmarek,F. Karsch, E. Laermann, C. Schmidt, and L. Scorzato,Phys. Rev. D , 074507 (2002).[15] D. Sexty, Phys. Lett. B , 108 (2014).[16] Y. Nambu and G. Jona-Lasinio, Phys. Rev. , 345(1961).[17] Y. Nambu and G. Jona-Lasinio, Phys. Rev. , 246(1961).[18] K. Fukushima, Physics Letters B , 277 (2004).[19] C. Ratti, M. A. Thaler, and W. Weise, Phys. Rev. D ,014019 (2006).[20] C. Sasaki, B. Friman, and K. Redlich, Phys. Rev. D ,074013 (2007).[21] S. Roessner, C. Ratti, and W. Weise, Phys. Rev. D ,034007 (2007).[22] K. Fukushima, Phys. Rev. D , 114028 (2008).[23] C. Sasaki, B. Friman, and K. Redlich, Phys. Rev. Lett. , 232301 (2007).[24] C. Sasaki, B. Friman, and K. Redlich, Phys. Rev. D ,054026 (2007).[25] B.-J. Schaefer, J. M. Pawlowski, and J. Wambach, Phys.Rev. D , 074023 (2007).[26] B.-J. Schaefer, M. Wagner, and J. Wambach, Phys. Rev.D , 074013 (2010).[27] A. J. Mizher, M. N. Chernodub, and E. S. Fraga, Phys.Rev. D , 105016 (2010).[28] T. K. Herbst, J. M. Pawlowski, and B.-J. Schaefer, Phys.Lett. B , 58 (2011).[29] V. Skokov, B. Stokic, B. Friman, and K. Redlich, Phys.Rev. C , 015206 (2010).[30] V. Skokov, B. Friman, E. Nakano, K. Redlich, and B. J.Schaefer, Phys. Rev. D , 034029 (2010).[31] V. Skokov, B. Friman, and K. Redlich, Phys. Rev. C ,054904 (2011).[32] M. Mitter and B.-J. Schaefer, Phys. Rev. D , 054027(2014). [33] T. K. Herbst, M. Mitter, J. M. Pawlowski, B.-J. Schaefer,and R. Stiele, Phys. Lett. B , 248 (2014).[34] M. Mitter, J. M. Pawlowski, and N. Strodthoff, Phys.Rev. D , 054035 (2015).[35] J. Braun, L. Fister, J. M. Pawlowski, and F. Rennecke,Phys. Rev. D , 034016 (2016).[36] C. S. Fischer, J. Luecker, and C. A. Welzbacher, Phys.Rev. D , 034022 (2014).[37] M. M. Aggarwal et al. (STAR), (2010), arXiv:1007.2613[nucl-ex].[38] P. Braun-Munzinger, K. Redlich, and J. Stachel,arXiv:nucl-th/0304013 [nucl-th].[39] A. Andronic, P. Braun-Munzinger, and J. Stachel, Nucl.Phys. A , 167 (2006).[40] A. Andronic, P. Braun-Munzinger, and J. Stachel, Phys.Lett. B , 142 (2009).[41] J. Cleymans, H. Oeschler, K. Redlich, and S. Wheaton,Phys. Rev. C , 034905 (2006).[42] M. A. Stephanov, K. Rajagopal, and E. V. Shuryak,Phys. Rev. Lett. , 4816 (1998).[43] M. A. Stephanov, K. Rajagopal, and E. V. Shuryak,Phys. Rev. D , 114028 (1999).[44] M. Asakawa, U. W. Heinz, and B. Muller, Phys. Rev.Lett. , 2072 (2000).[45] S. Jeon and V. Koch, Phys. Rev. Lett. , 2076 (2000).[46] B. Friman, F. Karsch, K. Redlich, and V. Skokov, TheEuropean Physical Journal C , 1694 (2011).[47] S. Ejiri, C. R. Allton, M. Doring, S. J. Hands, O. Kacz-marek, F. Karsch, E. Laermann, and K. Redlich, Nucl.Phys. Proc. Suppl. , 505 (2005).[48] C. R. Allton, M. Doring, S. Ejiri, S. J. Hands, O. Kacz-marek, F. Karsch, E. Laermann, and K. Redlich, Phys.Rev. D , 054508 (2005).[49] S. Ejiri, F. Karsch, and K. Redlich, Phys. Lett. B ,275 (2006).[50] F. Karsch, S. Ejiri, and K. Redlich, Nucl. Phys. A ,619 (2006).[51] F. Karsch and K. Redlich, Phys. Lett. B , 136 (2011).[52] L. Adamczyk et al. (STAR Collaboration), Phys. Rev.Lett. , 032302 (2014).[53] X. Luo (STAR), PoS CPOD2014 , 019 (2015).[54] X. Luo, Nuclear Physics A , 75 (2016).[55] A. Rustamov, Talk at Quark Matter 2017 .[56] C. Wetterich, Phys. Lett. B , 90 (1993).[57] T. R. Morris, Int. J. Mod. Phys. A , 2411 (1994).[58] J. Berges, N. Tetradis, and C. Wetterich, Phys. Rept. , 223 (2002).[59] J. Polonyi, Central Eur. J. Phys. , 1 (2003).[60] A. Bazavov et al. , Phys. Rev. Lett. , 192302 (2012).[61] P. M. Lo, B. Friman, O. Kaczmarek, K. Redlich, andC. Sasaki, Phys. Rev. D , 074502 (2013).[62] F. Rennecke, Phys. Rev. D , 076012 (2015).[63] J. Eser, M. Grahl, and D. H. Rischke, Phys. Rev. D ,096008 (2015).[64] D. F. Litim, Phys. Rev. D , 105007 (2001).[65] D. F. Litim, Phys. Lett. B , 92 (2000).[66] K. Kamikado, N. Strodthoff, L. von Smekal, andJ. Wambach, Phys. Lett. B , 1044 (2013).[67] R.-A. Tripolt, N. Strodthoff, L. von Smekal, andJ. Wambach, Phys. Rev. D , 034010 (2014).[68] J. P. Boyd, Chebyshev and Fourier Spectral Methods, Sec- ond Edition (DOVER Publications, Inc., 2001).[69] J. Borchardt and B. Knorr, Phys. Rev. D , 105011(2015).[70] J. Borchardt and B. Knorr, Phys. Rev. D , 025027(2016).[71] G. Almasi, R. Pisarski, and V. Skokov, (2016),arXiv:1612.04416 [hep-ph].[72] In the following, ω represents the expectation value ofthe Euclidean zero-component field, which differs fromthe Minkowski value by a factor of i .[73] O. Kaczmarek, F. Karsch, E. Laermann, C. Miao,S. Mukherjee, et al. , Phys. Rev. D , 014504 (2011).[74] H.-T. Ding, A. Bazavov, F. Karsch, Y. Maezawa,S. Mukherjee, et al. , Proc. Sci. LATTICE2013 , 157(2014).[75] H.-T. Ding and P. Hegde (Bielefeld-BNL-CCNU), Proc.Sci. LATTICE2015 , 161 (2016).[76] Y. Hatta and T. Ikeda, Phys. Rev. D , 014028 (2003).[77] A. Bazavov et al. , Phys. Rev. D , 054504 (2017).[78] M. Kitazawa, T. Koide, T. Kunihiro, and Y. Nemoto,Prog. Theor. Phys. , 929 (2002).[79] Since the CEP is a singular point, the correspondingvalue of χ B depends on how the limit is taken. Owingto the curvature of the phase boundary, the CEP is ap-proached slightly from below, which implies that χ B di-verges to plus infinity. [80] F. Karsch, J. Phys. Conf. Ser. , 012015 (2017).[81] A. Bzdak and V. Koch, Phys. Rev. C , 044904 (2012).[82] A. Bzdak, V. Koch, and V. Skokov, Phys. Rev. C ,014901 (2013).[83] M. Kitazawa and M. Asakawa, Phys. Rev. C , 021901(2012).[84] M. Kitazawa and M. Asakawa, Phys. Rev. C , 024904(2012).[85] G. A. Alm´asi, B. Friman, and K. Redlich, NuclearPhysics A , 356 (2016).[86] W.-j. Fu and J. M. Pawlowski, Phys. Rev. D , 091501(2016).[87] R. Esha, Talk at Quark Matter 2017 .[88] V. Skokov, B. Friman, and K. Redlich, Phys. Rev. C ,034911 (2013).[89] P. Braun-Munzinger, A. Rustamov, and J. Stachel, Nu-clear Physics A , 114 (2017).[90] L. F. Palhares, E. S. Fraga, and T. Kodama, J. Phys. G , 085101 (2011).[91] E. S. Fraga, L. F. Palhares, and P. Sorensen, Phys. Rev.C , 011903 (2011).[92] L. F. Palhares and E. S. Fraga, Phys. Atom. Nucl. ,906 (2012).[93] M. Hippert, E. S. Fraga, and E. M. Santos, Phys. Rev.D , 014029 (2016).[94] T. K. Herbst, J. M. Pawlowski, and B.-J. Schaefer, Phys.Rev. D88