Baryon number fluctuations in chiral effective models and their phenomenological implications

TL;DR

This paper investigates net-baryon-number fluctuations near the chiral phase transition using a QCD-inspired model, comparing results with experimental data and exploring effects of model assumptions and nonequilibrium dynamics.

Contribution

It applies the Functional Renormalization Group to a Polyakov-loop extended Quark-Meson model to analyze baryon fluctuations and compares theoretical predictions with RHIC data.

Findings

Model results agree with experimental data at higher energies.

Discrepancy observed in the fourth cumulant below 19.6 GeV.

Vector interactions influence the location of the critical endpoint.

Abstract

We study the critical properties of net-baryon-number fluctuations at the chiral restoration transition in a medium at finite temperature and net baryon density. The chiral dynamics of quantum chromodynamics (QCD) is modeled by the Polykov-loop extended Quark-Meson Lagrangian, that includes the coupling of quarks to vector meson and temporal gauge fields. The Functional Renormalization 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 $\chi_B^n$, for ${1\leq n\leq 6}$, near the phase boundary. The results are presented in the context of the recent experimental data of the STAR Collaboration on fluctuations of the net proton number in heavy ion collisions at RHIC. We show that the model results for the energy dependence of the cumulant ratios are in good overall agreement with the data, with one exception. At center-of-mass energies below ${19.6\;\mathrm{GeV}}$, we find that the measured fourth-order cumulant deviates considerably from the model results, which incorporate the expected $O(4)$ and $Z(2)$ criticality. We assess the influence of model assumptions and in particular of repulsive vector-interactions, which are used to modify the location of the critical endpoint in the model, on the cumulants ratios. Finally, we discuss a possibility to test to what extent the fluctuations are affected by nonequilibrium dynamics by comparing certain ratios of cumulants.

Tables

Table 1. Table 1: Parameters of the Polyakov-loop potential introduced in Eqs. ( 30 ) and ( 32 ).

$a_1$	$a_2$	$a_3$	$a_4$	$a_5$
-44.14	151.4	-90.0677	2.77173	3.56403
$b_1$	$b_2$	$b_3$	$b_4$
-0.32665	-82.9823	3.0	5.85559
$c_1$	$c_2$	$c_3$	$c_4$	$c_5$
-50.7961	114.038	-89.4596	3.08718	6.72812
$d_1$	$d_2$	$d_3$	$d_4$	$d_5$
27.0885	-56.0859	71.2225	2.9715	6.61433

Table 2. Table 2: Different combinations of the initial conditions and the Polyakov-loop potentials, used in the PQM model.

SET	Initial conditions	Polyakov-loop potential
A	(35)	(30)
B	(35)	(28)
C	(36)	(30)
D	(36)	(28)

Equations

$$\Phi$$

$$\bar{\Phi }$$

$$\mathcal{L}=$$

$$+\frac{1}{2}{\left({\partial }_{\mu }\sigma \right)}^2+\frac{1}{2}{\left({\partial }_{\mu }\vec{\pi }\right)}^2-U_m(\sigma ,\vec{\pi })$$

$$-\mathcal{U}(\Phi ,\bar{\Phi };T)-\frac{1}{2}{m}_{\omega }^2{\omega }^2+\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu },$$

$$U_m(\sigma ,\vec{\pi })=\frac{\lambda }{4}({\sigma }^2+{\vec{\pi }}^2{)}^2+\frac{m^2}{2}({\sigma }^2+{\vec{\pi }}^2)-H\sigma ,$$

$${\partial }_k{\Gamma }_k[\varphi ]=\frac{1}{2}\mathrm{STr}\left[{\left({\Gamma }_k^{(2)}[\varphi ]+R_k\right)}^{-1}{\partial }_k{R}_k\right],$$

$$\Omega (T,{\mu }_q)=-P(T,{\mu }_q)={\Omega }_{k=0}(T,{\mu }_q,{\sigma }^{\ast },{\Phi }^{\ast },{\bar{\Phi }}^{\ast })-H{\sigma }^{\ast },$$

$${\partial }_k$$

$$=$$

$$\displaystyle-\frac{4N_{c}N_{f}}{E_{q}}\big{(}1-N(T,\mu_{q};\sigma,\Phi,\bar{\Phi})-N(T,-\mu_{q};\sigma,\bar{\Phi},\Phi)\big{)}\bigg{\}},$$

$$N$$

$$=\frac{1+2\bar{\Phi }{e}^{(E_q-{\mu }_q)/T}+\Phi e^{2(E_q-{\mu }_q)/T}}{1+3\bar{\Phi }{e}^{(E_q-{\mu }_q)/T}+3\Phi e^{2(E_q-{\mu }_q)/T}+e^{3(E_q-{\mu }_q)/T}}.$$

$${\Omega }_{k=\Lambda }={U_m(\sigma ,0)|}_{H=0}+\mathcal{U}(\Phi ,\bar{\Phi };T).$$

$${\frac{\partial {\Omega }_{k=0}}{\partial \sigma }|}_{SP}=H,{\frac{\partial {\Omega }_{k=0}}{\partial \Phi }|}_{SP}=0,{\frac{\partial {\Omega }_{k=0}}{\partial \bar{\Phi }}|}_{SP}=0.$$

$$\frac{\partial {\Omega }_k}{\partial \Phi },$$

$$\frac{\partial {\Omega }_k^2}{\partial \mu \partial \Phi },$$

$${\Omega }_{k=\Lambda }={U_m(\sigma ,0)|}_{H=0}+\mathcal{U}(\Phi ,\bar{\Phi };T)-\frac{1}{2}{m}_{\omega }^2{\omega }^2.$$

$${\frac{\partial {\Omega }_{k=0}}{\partial \omega }|}_{SP}=0,⟹g_{\omega }\omega =-G_{\omega }\frac{\partial \tilde{\Omega }}{\partial \nu }=G_{\omega }\tilde{n}(T,\nu ),$$

$$\mu =\nu +G_{\omega }\tilde{n}(T,\nu ).$$

$$\Omega (T,\mu )=-P(T,\mu )=\tilde{\Omega }(T,\nu )-\frac{(\mu -\nu {)}^2}{2G_{\omega }}$$

$$3^n{\chi }_B^n={\chi }_q^n=\frac{1}{T^{4-n}}\frac{{\partial }^n\Omega (T,{\mu }_q)}{\partial {\mu }_q^n}.$$

$$\displaystyle\chi_{B}^{2}=\frac{1}{9T^{2}}\bigg{(}$$

$$\displaystyle-\sum_{i,j}\frac{\partial^{2}\Omega_{k=0}}{\partial\mu_{q}\partial\varphi_{i}}\left(\frac{\partial^{2}\Omega_{k=0}}{\partial\varphi\partial\varphi}\right)^{-1}_{ij}\frac{\partial^{2}\Omega_{k=0}}{\partial\mu_{q}\partial\varphi_{j}}\bigg{)},$$

$$\frac{d\nu }{d\mu }=\frac{1}{1+9T^2{G}_{\omega }{\tilde{\chi }}_B^2(T,\nu )},$$

$$\begin{array}{rl}{\chi }_B^1& ={\tilde{\chi }}_B^1,\\ {\chi }_B^2& =\frac{{\tilde{\chi }}_B^2}{1+9T^2{G}_{\omega }{\tilde{\chi }}_B^2(T,\nu )},\\ {\chi }_B^3& =\frac{{\tilde{\chi }}_B^3}{{\left(1+9T^2{G}_{\omega }{\tilde{\chi }}_B^2\right)}^3},\\ {\chi }_B^4& =\frac{{\tilde{\chi }}_B^4}{{\left(1+9T^2{G}_{\omega }{\tilde{\chi }}_B^2\right)}^4}-\frac{27T^2{G}_{\omega }{\left({\tilde{\chi }}_B^3\right)}^2}{{\left(1+9T^2{G}_{\omega }{\tilde{\chi }}_B^2\right)}^5}.\end{array}$$

$${\chi }_B^{n,m}$$

$$\kappa {\sigma }^2$$

$${\chi }_B^{2m-1,2n-1}={\chi }_B^{2m,2n}$$

$$\frac{d}{dk}{\Omega }_k=3F_b\left(m_{\pi }^2\right)+F_b\left(m_{\sigma }^2\right)+F_f\left(\mu ,\Phi ,\bar{\Phi }\right),$$

$$F_b(m^2)$$

$$F_f(\mu ,\Phi ,\bar{\Phi })$$

$$m_{\pi }^2=\frac{1}{\sigma }\frac{\partial {\Omega }_k}{\partial \sigma },m_{\sigma }^2=\frac{{\partial }^2{\Omega }_k}{\partial {\sigma }^2}.$$

$$\begin{array}{rl}\frac{d}{dk}\frac{\partial {\Omega }_k}{\partial {\lambda }_i}& =\frac{\partial }{\partial {\lambda }_i}\frac{d{\Omega }_k}{dk}=3F_b^{\prime }\left(m_{\pi }^2\right)\frac{1}{\sigma }\frac{{\partial }^2{\Omega }_k}{\partial {\lambda }_i\partial \sigma }+F_b^{\prime }\left(m_{\sigma }^2\right)\frac{{\partial }^3{\Omega }_k}{\partial {\lambda }_i{\partial }^2\sigma }+\frac{\partial }{\partial {\lambda }_i}{F}_f(\mu ,\Phi ,\bar{\Phi }),\\ \frac{d}{dk}\frac{{\partial }^2{\Omega }_k}{\partial {\lambda }_i\partial {\lambda }_j}& =\frac{{\partial }^2}{\partial {\lambda }_i\partial {\lambda }_j}\frac{d{\Omega }_k}{dk}=3F_b^{\prime \prime }\left(m_{\pi }^2\right)\frac{1}{{\sigma }^2}\frac{{\partial }^2{\Omega }_k}{\partial {\lambda }_i\partial \sigma }\frac{{\partial }^2{\Omega }_k}{\partial {\lambda }_j\partial \sigma }+3F_b^{\prime }\left(m_{\pi }^2\right)\frac{1}{\sigma }\frac{{\partial }^3{\Omega }_k}{\partial {\lambda }_i\partial {\lambda }_j\partial \sigma }\\ & +F_b^{\prime \prime }\left(m_{\sigma }^2\right)\frac{{\partial }^3{\Omega }_k}{\partial {\lambda }_i{\partial }^2\sigma }\frac{{\partial }^3{\Omega }_k}{\partial {\lambda }_j{\partial }^2\sigma }+F_b^{\prime }\left(m_{\sigma }^2\right)\frac{{\partial }^4{\Omega }_k}{\partial {\lambda }_i\partial {\lambda }_j{\partial }^2\sigma }+\frac{{\partial }^2}{\partial {\lambda }_i\partial {\lambda }_j}{F}_f(\mu ,\Phi ,\bar{\Phi }),\end{array}$$

$$\frac{\mathcal{U}(\Phi ,\bar{\Phi };T)}{T^4}=-\frac{b_2(T)}{2}\left(\Phi \bar{\Phi }\right)-\frac{b_3}{6}\left({\Phi }^3+{\bar{\Phi }}^3\right)+\frac{b_4}{4}{\left(\Phi \bar{\Phi }\right)}^2,$$

Full text

I Introduction

It is well established within lattice QCD (LQCD), that at vanishing and small chemical potential, $\mu_{B}$, strongly interacting matter undergoes a smooth crossover from the hadronic phase to a quark-gluon plasma Aoki et al. (2006, 2009); Borsanyi et al. (2010); Bazavov et al. (2012a). Moreover, arguments are given, that in the limit of massless $u$ and $d$ quarks, the crossover transition becomes a genuine second-order chiral phase transition, belonging to the $O(4)$ universality class Pisarski and Wilczek (1984); Ejiri et al. (2009). The nature of this transition at higher net baryon densities is still not settled by first principle LQCD studies, owing to the sign problem. However, in effective models of QCD, it is found that at sufficiently large $\mu_{B}$, the system can exhibit a first order chiral phase transition. The endpoint of this conjectured transition line in the $(T,\mu_{B})$- plane, is the chiral critical endpoint (CEP) Asakawa and Yazaki (1989); Barducci et al. (1989); Wilczek (1992); Halasz et al. (1998); Ding et al. (2015). At the CEP, the system exhibits the 2nd order phase transition, which belongs to the $Z(2)$ universality class Schaefer and Wambach (2007).

Various methods have been developed to circumvent the sign problem and perform LQCD calculations at nonvanishing $\mu_{B}$ Fodor and Katz (2002); Allton et al. (2002); Sexty (2014). However, so far these methods are restricted to rather small values of chemical potential, ${\mu_{B}\ll 3T_{c}}$, where ${T_{c}=\;154(9)\mathrm{MeV}}$ is the chiral crossover temperature at vanishing net baryon density Borsanyi et al. (2010); Bazavov et al. (2012a).

Due to the restriction of present LQCD calculations to small net baryon densities, effective models that belong to the same universality class as QCD, e.g., the Polyakov-loop extended Nambu-Jona-Lasinio Nambu and Jona-Lasinio (1961a, b); Fukushima (2004); Ratti et al. (2006); Sasaki et al. (2007a); Roessner et al. (2007); Fukushima (2008); Sasaki et al. (2007b, c) and Quark-Meson (PQM) models Schaefer et al. (2007, 2010); Mizher et al. (2010); Herbst et al. (2011); Skokov et al. (2010a, b, 2011); Mitter and Schaefer (2014); Herbst et al. (2014), have been employed to study the chiral phase transition for a broad range of thermal parameters. Moreover, there are ongoing studies using functional methods that focus on systematic improvements of such effective models, to achieve a more quantitative matching with LQCD, and eventually to perform reliable calcuations of the properties of QCD matter at nonzero net baryon density Mitter et al. (2015); Braun et al. (2016); Fischer et al. (2014).

One of the strategic goals of current experimental and theoretical studies of chiral symmetry restoration in QCD is to unravel the phase diagram of strongly interacting matter and to clarify whether a chiral CEP exists. A dedicated research program at RHIC, the beam energy scan (BES), has been established to explore these issues in collisions of heavy ions at relativistic energies Aggarwal et al. (2010). By varying the beam energy at RHIC, the properties of strongly interacting matter in a broad range of net baryon densities, correponding to a wide range in baryon chemical potential, ${20\;\mathrm{MeV}<\mu_{B}<500\;\mathrm{MeV}}$ Braun-Munzinger et al. ; Andronic et al. (2006, 2009); Cleymans et al. (2006), can be studied. In particular, fluctuations of conserved charges are considered as relevant probes of the phase diagram Stephanov et al. (1998, 1999); Asakawa et al. (2000); Jeon and Koch (2000); Friman et al. (2011); Ejiri et al. (2005); Allton et al. (2005); Ejiri et al. (2006); Karsch et al. (2006); Karsch and Redlich (2011); Friman et al. (2011); Sasaki et al. (2007b, c). These are experimentally accessible and reflect the criticality of the chiral transition. Fluctuations of the net baryon number are particulary interesting, owing to a direct connection to critical scaling near the chiral phase boundary Friman et al. (2011); Karsch and Redlich (2011).

First data on net-proton-number fluctuations, which are used as a proxy for fluctuations of the net baryon number, in heavy-ion collisions have been obtained by the STAR Collaboration at RHIC energies Adamczyk et al. (2014); Luo (2015, 2016). The STAR data on the variance, skewness and kurtosis of net proton number, are intriguing and have stimulated lively discussions on their physics origin and interpretation. Recently, first results on the mean and variance of the net-proton distribution were obtained in heavy ion collisions at the LHC energy by the ALICE Collaboration Rustamov .

In the following, we focus on the properties and systematics of the cumulants of net-baryon-number fluctuations near the chiral phase boundary. Our studies are performed based on the Polykov-loop extended Quark-Meson Lagrangian which includes couplings of quarks to vector and temporal gauge fields. To account for the ${O(4)}$ and ${Z(2)}$ critical fluctuations near the phase boundary, we employ the Functional Renormalisation Group (FRG) Wetterich (1993); Morris (1994); Berges et al. (2002); Polonyi (2003). We present contour plots of ratios of cumulants, involving the mean, the variance, skewness and kurtosis in the ${(T,\mu_{B})}$- plane. These ratios obtained on the phase boundary and on a freeze-out line determined by fitting the skewness ratio, following Ref. Bazavov et al. (2012b), are confronted with the corresponding experimental values of the STAR Collaboration. We explore the influence of the CEP and the repulsive interactions on the fluctuation observables. The calculations are performed with two different initial conditions and two Polyakov-loop potentials, to assess the model dependence of the results.

We find that the energy dependence of ratios of low-order cumulants, $\chi_{B}^{n}$ with ${n=1,2,3}$, are in good agreement with the data, whereas for the ratio involving the kurtosis, the model results differ from the data at energies below $\sqrt{s}\simeq 20$ GeV, corresponding roughly to the top SPS energy. At these low energies, the latter ratio increases strongly, while the model results, which embody $O(4)$ as well as $Z(2)$ criticality, differ substantially from the data. Finally, we discuss possible caveats, which could undermine our conclusions and assess the uncertainties of the model.

The paper is organized as follows. In Sec. II we formulate the model and the FRG method, which is employed to compute cumulants of net baryon number. In Sec. III we present results on the characteristics of susceptibility ratios as functions of temperature and chemical potential and study, along the phase boundary, the dependence of these ratios on the vector interaction. Moreover, in this section we confront the model results with the STAR data on net-proton-number fluctuations. Sec. IV contains summary and conclusions.

II The Polyakov–quark-meson model

The PQM model is a low energy effective approximation to QCD formulated in terms of the light quark $q=(u,d)$ as well as scalar and the pseudoscalar meson $\phi=(\sigma,\vec{\pi})$ fields. The quarks are coupled to the background Euclidean gluon field $A_{\mu}$, with vanishing spatial components, which is linked to the Polyakov loop

[TABLE]

Moreover, we include a coupling of the quarks to a massive vector field $\omega$. The resulting Lagrangian of the model reads

[TABLE]

where $D_{\mu}=\partial_{\mu}-iA_{\mu}$, and $F^{\mu\nu}=\partial^{\mu}\omega^{\nu}-\partial^{\nu}\omega^{\mu}$. The parameters of the mesonic potential

[TABLE]

are tuned to vacuum properties of the $\sigma$ and $\vec{\pi}$ mesons (see Appendix B). We use the Polyakov-loop potential $\mathcal{U}(\Phi,\bar{\Phi};T)$ determined in Lo et al. (2013), by fitting quenched lattice QCD results for the equation of state and the Polyakov-loop susceptibilities. The parametrization of ${\mathcal{U}(\Phi,\bar{\Phi};T)}$ is given in Appendix B.

We compute the thermodynamic properties of this model including fluctuations of the scalar and pseudoscalar meson fields within the framework of the FRG method. The $\omega$ meson, on the other hand, is treated on the mean-field level. The MF treatment of the vector field is justified by recent FRG results obtained at vanishing chemical potential indicating that the vector meson mass remains above the cutoff during the FRG flow and that the temperature dependence of the screening mass is very weak Rennecke (2015); Eser et al. (2015). Therefore, it is expected that fluctuations of the vector fields decouple from the flow. The Polyakov loop is also treated on the mean-field level. The Polyakov-loop variables and the vector field are tuned such that at the end of the calculation a stationary point of the thermodynamic potential is reached.

In the FRG framework, the effective average action $\Gamma_{k}$, which interpolates between the classical and the full quantum action, is obtained by solving the renomalization group flow equation Wetterich (1993)

[TABLE]

where $\phi$ denotes the quantum fields considered, $\mathrm{STr}$ is a trace over the fields, over momentum and over all internal indices. It also adds fermion contribution to the boson contribution with correct sign and prefactor. The regulator function $R_{k}$ suppresses fluctuations at momenta below $k$. Thus, effects of fluctuations of quantum fields are included shell by shell in momentum space, starting from a UV cutoff scale $\Lambda$. We employ the optimized regulator introduced by Litim Litim (2001, 2000), which yields an algebraic expression for the right hand side of the flow equation, Eq. (5).

II.1 Thermodynamics at vanishing vector coupling

For vanishing vector coupling, the vector fields obviously decouple completely. Hence, in such models the thermodynamics is controlled by the quark, as well as, the scalar and pseudoscalar fields Herbst et al. (2011); Skokov et al. (2010a); Herbst et al. (2014); Kamikado et al. (2013); Tripolt et al. (2014).

Using $\left.\Omega_{k}=\Gamma_{k}/V\right|_{H=0}$ we obtain the thermodynamic potential density

[TABLE]

where $\mu_{q}=\mu_{B}/3$ denotes the quark chemical potential, the asterisk denotes the expectation value of the condensates and $\Omega_{k=0}$ is obtained by solving the flow equation,

[TABLE]

with $E_{q}=\sqrt{k^{2}+g_{s}^{2}\sigma^{2}}$, $E_{\sigma}=\sqrt{k^{2}+\Omega_{k}^{\prime\prime}}$ and $E_{\pi}=\sqrt{k^{2}+\Omega_{k}^{\prime}/\sigma}$. Here the prime denotes a derivative with respect to $\sigma$ field, $g_{s}$ is the Yukawa coupling, while $N_{c}=3$ is the number of colors and $N_{f}=2$ the number of flavors. The Polyakov-loop modified quark occupation numbers are given by

[TABLE]

The flow equation (II.1) is solved numerically starting from the initial conditions specified at the momentum scale $k=\Lambda$

[TABLE]

The thermodynamic potential $\Omega(T,\mu)$, shown in Eq. (6), is obtained by evaluating the fully evolved $\Omega_{k=0}$ at the stationary point, $SP=(\sigma^{*},\Phi^{*},\bar{\Phi}^{*})$, determined by

[TABLE]

Any thermodynamic quantity can be obtained by taking the appropriate derivatives of the thermodynamic potential (6). However, quantities that require higher order derivatives are difficult to extract by numerical differentiation, owing to limited numerical precision. We find that improved numerical stability is obtained when the derivatives up to second order,

[TABLE]

are computed directly using flow equations, summarized in Appendix A.

We solve the coupled flow equations numerically by discretizing them on a grid of the order parameter $\sigma$. In order to recover the Stefan-Boltzmann limit, it is necessary to include the effect of fermion fluctuations above the UV cutoff $\Lambda$. This is done by adding the contribution of the thermal Polyakov-loop suppressed massless quark loop, integrated from the UV cutoff to infinity Skokov et al. (2010a). The initial conditions for the flow equations are fixed by vacuum criteria, as discussed in Appendix B.

It should be noted that we also performed the RG calculations by expanding the potential in terms of Chebyshev polynomials Boyd (2001), instead of working on a uniformly discretized grid. This method is illustrated in Refs. Borchardt and Knorr (2015, 2016); Almasi et al. (2016). Using the same initial conditions it yields the same results for all thermodynamic quantities, but up to a higher numerical precision. This method was required to extract cumulants of high order.

II.2 Thermodynamics at nonvanishing vector coupling

For nonvanishing vector coupling, the $\omega$ field is coupled to the chiral sector through the fermionic part of the flow. The temporal component of the $\omega$ field gains a nonvanishing expectation value at nonzero density, while the spatial components vanish in a system with zero net quark current.

An inspection of the Lagrangian (II), reveals that a nonvanishing $\omega$ field effectively acts as a shift of the quark chemical potential. Thus, it is convenient to introduce an effective quark chemical potential, $\nu=\mu-g_{\omega}\omega$ 111In the following, $\omega$ represents the expectation value of the Euclidean zero-component field, which differs from the Minkowski value by a factor of $i$, and continue using the formalism for vanishing vector coupling, albeit with modified initial conditions for the flow equation

[TABLE]

Since, for a given value of $\omega$, this amounts to constant shift of the thermodynamic potential, it will not modify the RG flow. Thus, the only effect of a nonzero $\omega$ field on the flow is the replacement of the chemical potential $\mu$ by the previously defined effective chemical potential $\nu$.

The expectation value of $\omega$ is obtained by extremizing $\Omega_{k=0}$. Since the flow for vanishing and nonvanishing vector coupling can be related to each other, one can express all results in terms of quantities computed at vanishing vector coupling. We denote the thermodynamic potential density and net baryon density at vanishing vector coupling by $\tilde{\Omega}(T,\nu)$ and $\tilde{n}(T,\nu)$, respectively. The expectation value of the $\omega$ field is then determined by

[TABLE]

where $G_{\omega}=g_{\omega}^{2}/m_{\omega}^{2}$. This yields the relation between the real and the effective chemical potentials

[TABLE]

Moreover, the corresponding thermodynamic potential for nonzero vector coupling is given by

[TABLE]

We note, that for a first order phase transition, this procedure cannot be used to identify the critical value of thermal parameters, owing to a flattening of the potential, as a function of the order parameter.

II.3 Net-baryon-number cumulants

Fluctuations of conserved charges, in particular of the net baryon number, are valuable probes of critical phenomena and can be used to identify the location of the phase boundary of chiral symmetry restoration. The fluctuations are characterized by cumulants, that can be obtained directly from the partition function. The $n^{\mathrm{th}}$-order baryon $\chi_{B}^{n}$, and quark $\chi_{q}^{n}$ number cumulants are

[TABLE]

At vanishing vector coupling, the first order cumulant $\chi_{B}^{1}$, the net baryon density, can be obtained directly from the flow of $\partial\Omega_{k}/\partial\mu$, while the second order cumulant $\chi_{B}^{2}$ is given by

[TABLE]

where $\varphi=(\sigma,\Phi,\bar{\Phi})$. The derivatives on the right hand side of Eq. (II.3) are determined by the solution of flow equations introduced in Appendix A. Higher order cumulants are then computed using numerical differentiation of $\chi_{B}^{2}$.

For nonvanishing vector coupling, the baryon number cumulants can be expressed in terms of the corresponding cumulants $\tilde{\chi}_{B}^{n}$ for vanishing vector coupling. This is done using Eqs. (14) and (15). Based on Eq. (14) one can write

[TABLE]

which using Eq. (15), leads to the following relations

[TABLE]

Here the cumulants $\tilde{\chi}_{B}^{k}$ are functions of temperature and the reduced chemical potential $\nu$. Using these relations, we can compute the influence of the vector interaction on the baryon number cumulants, using the RG flow obtained for vanishing vector coupling, $G_{\omega}=0$.

III Net-baryon-number cumulants and the phase boundary

In the previous section we have introduced the method employed to compute cumulants of net baryon number in an effective chiral model in the presence of vector interactions. In the following, we discuss the properties of these cumulants near, and at the chiral phase boundary, at finite temperature and chemical potential.

The chiral Lagrangian introduced above, shares the chiral critical properties with QCD. In particular, at moderate values of the chemical potential, the PQM model exhibits a chiral transition belonging to the $O(4)$ universality class Pisarski and Wilczek (1984); Schaefer and Wambach (2007). For larger values of $\mu$, it reveals a $Z(2)$ critical endpoint, followed by a first order phase transition Schaefer and Wambach (2007). Consequently, the PQM model embodies the generic phase structure expected for QCD, with the universal $O(4)$ and $Z(2)$ criticality encoded in the scaling functions. Furthermore, due to the coupling of the quarks to the background gluon fields, the PQM model incorporates ”statistical confinement”, i.e., the suppression of quark and diquark degrees of freedom in the low temperature, chirally broken phase Fukushima (2004); Skokov et al. (2010b). Consequently, by studying fluctuations of conserved charges in the PQM model, one can explore the influence of chiral symmetry restoration and of ”statistical confinement” on the cumulants in different sections of the chiral phase boundary. This study is of particular interest in the context of heavy-ion collisions, where cumulants of conserved charges are expected to provide a characteristic signature for the QCD phase boundary and for the conjectured critical endpoint Stephanov et al. (1998, 1999); Asakawa et al. (2000); Jeon and Koch (2000); Friman et al. (2011); Ejiri et al. (2006); Karsch and Redlich (2011).

III.1 Criticality of net-baryon-number cummulants in the $(T,\mu)$ plane

Generalized susceptibilities of conserved charges, which are given by higher-order derivatives of the thermodynamic pressure with respect to the corresponding chemical potentials, may exhibit a nonmonotonic dependence on the thermodynamic parameters. This is particularly the case in the vicinity of phase boundary and the CEP. In the critical region of the chiral transition, the strength of the fluctuations and the sign of the susceptibilities are by and large determined by the singular part of the free energy, which is encoded in the universal scaling functions, common to QCD and the PQM model. Thus, generic structures of the susceptibilities and relations between them near the phase boundary, can also be studied in the PQM model. Of particular interest is the behavior of susceptibilities along the $O(4)$ crossover line, and their modification as the CEP is approached. In the model calculations, the position of the CEP depends on the strength of the vector interaction. Thus, by changing the vector coupling $g_{\omega}$, one can assess the dependence of critical fluctuations in different sectors of the phase boundary on the location of the CEP.

For nonvanishing light-quark masses, the chiral symmetry is explicitly broken, which implies that at finite temperatures and moderate values of the baryon chemical potential, the system exhibits a chiral crossover transition. Thus, at small $\mu_{q}$, the fluctuations of conserved charges remain finite. However, because QCD matter at the physical values of the $u$ and $d$ quark masses is within the critical region of the second order phase transition, the fluctuations are still influenced by $O(4)$ criticality Ejiri et al. (2009); Kaczmarek et al. (2011); Ding et al. (2014); Ding and Hegde (2016); Hatta and Ikeda (2003). At physical quark masses, a genuine phase transition in QCD and the associated singular behavior of fluctuations are only expected at the conjectured CEP and along the line of the first-order phase transition.

To explore the phase diagram with net-baryon-number fluctuations, we need to know the qualitative dependence of net-baryon-number cumulants on the temperature and baryon chemical potential. In Fig. 1 we show the $T$-dependence at $\mu_{q}=0$, and in Figs. 2 and  3, the contour plots in the $(T,\mu_{q})$-plane of ratios of net-baryon-number susceptibilities. We focus on the following ratios of net-baryon-number susceptibilities:

[TABLE]

where $M$ is the mean, $\sigma$ the variance, $S_{B}$ the skewness and $\kappa$ the kurtosis of the net-baryon-number distribution. The statistical descriptors are related to the susceptibilities through $M=VT^{3}\chi_{B}^{1}$, $\sigma^{2}=VT^{3}\chi_{B}^{2}$ etc..

At vanishing chemical potential, all odd susceptibilities of net baryon number vanish, owing to the baryon-antibaryon symmetry. In addition, in the $O(4)$ universality class, the second and fourth order cumulants remain finite at the phase transition temperature at $\mu_{q}=0$ even in the chiral limit, implying that only sixth and higher order susceptibilities diverge. Thus, for physical quark masses, only higher order cumulants, $\chi_{B}^{n}$ with $n>4$, can exhibit $O(4)$ criticality at $\mu_{q}=0$ Friman et al. (2011). A further consequence of the baryon-antibaryon symmetry is the equality of the ratios

[TABLE]

for any integer $m$ and $n\,\geq 1$ at $\mu_{q}=0$. For $\chi_{B}^{3,1}$ $\chi_{B}^{4,2}$, the equality at small $\mu_{q}$ can also be confirmed by a direct comparison of the right panel of Fig. 2 with the left panel of Fig. 3.

At finite net baryon density the singularity at the $O(4)$ critical line is stronger than at $\mu_{q}=0$. Thus, in this case the third-order cumulant and all higher-order ones, diverge at the $O(4)$ line. The second order cumulant $\chi_{B}^{2}$ remains finite, and only diverges at the tricritical point for vanishing quark masses, and at the CEP for nonzero quark masses.

In Fig. 1 we show results on the temperature dependence of several ratios $\chi_{B}^{n,m}$ of net-baryon-number susceptibilities, together with the variance of the chiral condensate, $\chi_{ch}$. The location of the maximum of the chiral susceptibility, $\chi_{ch}$, defines the pseudo-critical temperature, $T_{c}$.

As shown in Fig. 1, the kurtosis $\kappa\sigma^{2}$ exhibits a rapid drop by approximately an order of magnitude at the phase boundary. This strong reduction of the kurtosis is attributed to the deconfinement transition, where the degrees of freedom carrying baryon number change from baryons to quarks Karsch et al. (2006); Ejiri et al. (2006). A qualitative understandning of this issue is obtained in the Boltzmann approximation, where the baryon contribution to the thermodynamic pressure is of the form $P_{B}(T,\mu_{B})\simeq F(T)\cosh(B\mu_{B}/T)$. Here $B$ is the baryon number of the relevant degrees of freedom and $F(T)$ a function of the corresponding excitation spectrum. At low temperatures, the degrees of freedom are baryons, with $B=\pm 1$, while at high temperatures they are quarks with $B=\pm 1/3$. Consequently, the kurtosis ratio, $\kappa\sigma^{2}=\chi_{B}^{4}/\chi_{B}^{2}$, is approximately proportional to the square of the baryon number of the relevant degrees of freedom Karsch et al. (2006); Ejiri et al. (2006). This clarifies the cause for the rapid change of $\kappa\sigma^{2}$ at the crossover transition from $\simeq 1$ to $\simeq 1/9$, seen in Fig. 1. In the limit where the characteristic mass is small compared to the temperature, the kurtosis ratio is modified by quantum statistics, to $\kappa\sigma^{2}=(1/9)(6/\pi^{2})$. We note, that the kurtosis ratio is independent of the mass spectrum as well as of any kinematic cuts, as long as the Boltzmann limit remains valid.

The ratio $\chi_{B}^{6,2}$, also shown in Fig. 1, exhibits more structure at the chiral transition. The observed characteristic temperature dependence of this ratio, in particular with a region of negative values at $T\simeq T_{c}$, is a consequence of the residual chiral $O(4)$ criticality Friman et al. (2011). In the absence of chiral critical fluctuations, $\chi_{B}^{6,2}$ would show a similar behavior as $\kappa\sigma^{2}$, with a smooth reduction from unity at low temperatures towards zero above the deconfinemnt transition. The distinctive behavior of the $\chi_{B}^{6,2}$ and $\chi_{B}^{4,2}$ ratios was already obtained in the PQM model within the FRG approach Friman et al. (2011), and agrees qualitatively with LQCD results Bazavov et al. (2017). Thus, the PQM model correctly captures the physics of QCD related to deconfinement and to the critical dynamics at the chiral transition. The negative values of $\chi_{B}^{6,2}$ at the chiral crossover were proposed as a signature for partial restoration of chiral symmetry in heavy ion collisions Friman et al. (2011).

As shown in Fig. 1, in the presence of repulsive interactions, $\chi_{B}^{6,2}$ and the kurtosis ratios are modified. The vector interaction leads to a downward shift in temperature of the ratios relative to the phase boundary, as well as, a suppression of the sixth order susceptibility in the temperature range below $T_{c}$. Nevertheless, the qualitative form of the ratios is preserved. In particular, the characteristic structure, where the sixth order cumulant is negative in a range of temperatures near $T_{c}$ owing to $O(4)$ criticality, is not eliminated by the repulsive interaction.

In Fig. 2 we show contour plots of the ratios $\chi_{B}^{1,2}$ and $\chi_{B}^{3,1}$ in the $(T,\mu_{q})$-plane. As noted above, all odd cumulants vanish at $\mu_{q}=0$, owing to baryon-antibaryon symmetry. Consequently, $\chi_{B}^{1,2}|_{\mu_{q}=0}=0$ for any $T$, while the ratio $\chi_{B}^{3,1}|_{\mu_{q}=0}$ is nonvanishing. At low and high $T$, relative to $T_{c}$, this ratio is consistent with unity and $2/(3\pi^{2})$, respectively, as expected since $\chi_{B}^{3,1}|_{\mu_{q}=0}=\chi_{B}^{4,2}|_{\mu_{q}=0}$. As indicated in Figs. 1 and 2, the ratio $\chi_{B}^{3,1}$ decreases with temperature, and depends weakly on the chemical potential. Moreover, Fig. 2, shows that for $\mu_{q}<T$, $\chi_{B}^{1,2}$ increases with $\mu_{q}$, and is only weakly dependent on the temperature. Thus, the ratio $\chi_{B}^{3,1}$ can be used as a measure of the temperature, while $\chi_{B}^{1,2}$ provides a gauge of the chemical potential.

At small $\mu_{q}/T$, the properties of the first four susceptibilities, $\chi_{B}^{n}$ with $n=1,..,4$, and consequently their ratios near the chiral crossover are dominantly affected by the coupling of the quarks to the Polyakov loop, and the resulting statistical confinement. The critical chiral dynamics, i.e. the $O(4)$ and $Z(2)$ criticality at the chiral crossover transition and at the CEP, respectively, unfolds at larger $\mu_{q}/T$. Near the CEP, there is a strong variation of the cumulants with $T$ and $\mu_{q}$, which increases with the order of cumulants. As shown in Figs. 2 and 3, the qualitative behavior of the cumulant ratios on lines going to the CEP depends strongly on the direction from which the CEP is approached. This behavior is governed by the critical properties encoded in the $O(4)$ and $Z(2)$ universal scaling functions.

The influence of the CEP on the characteristics of the various cumulant ratios can be studied by varying the value of the vector coupling, thus changing the position of the CEP in the $(T,\mu_{q})$-plane. This is illustrated in Fig. 3 for the kurtosis ratio. A comparison of the results shown in the left and right contour plots shows that with increasing $g_{\omega}$, the curvature of the phase boundary is reduced and the position of the CEP is shifted to lower $T$ and larger $\mu_{q}$ Kitazawa et al. (2002). As is clearly seen, when comparing the left and right plots of Fig. 3, there is a corresponding shift of the contours of $\chi_{B}^{4,2}$. These results clearly illustrate the influence of the critical endpoint on the fluctuation observables. A shift of the CEP to larger net baryon density, suppresses the magnitude of the net-baryon-number susceptibilities at a given $T$ and $\mu_{q}$.

III.2 Net-baryon-cumulant ratios and freeze-out in heavy ion collisions

In heavy-ion collisions, the thermal fireball formed in the quark-gluon plasma phase undergoes expansion and passes through the QCD phase boundary at some point $(T_{c},\mu_{c})$ point, which depends on the collision energy, $\sqrt{s}$. Analysis of ratios of particle multiplicities indicate that at high beam energies (small values of $\mu_{q}/T$), the freeze-out occurs at, or just below the phase boundary. Thus, the beam energy dependence of net-baryon-number susceptibilities can provide insight into the structure of the QCD phase diagram and information on the existence and location of the CEP. Consequently, it is of phenomenological interest to compute the properties of fluctuations of conserved charges along the chiral phase boundary. Since there, the critical structure and the relations between different susceptibilities are governed by the universal scaling functions, the generic behavior of ratios of net-baryon-number susceptibilities can be explored also in model calculations.

In Fig. 4 we show $\chi_{B}^{1,2}$, as well as skewness and kurtosis ratios of net baryon number, obtained in the PQM model along the chiral phase boundary, with and without vector repulsion. In the former case, the position of the CEP is shifted to smaller $T$ and larger $\mu_{q}$.

The ratio $\chi_{B}^{1,2}$ exhibits a maximum along the phase boundary. For small $(\mu_{q}/T)_{c}$ this ratio is well approximated by $\tanh(3\mu_{q}/T)$, while after reaching a maximum, it decreases as the CEP is approached. The cumulants $\chi_{B}^{1}$ and $\chi_{B}^{2}$ remain finite along the $O(4)$ line. However, at the CEP the variance of the net-baryon-number fluctuations diverges. This is the cause for the observed decrease of $\chi_{B}^{1,2}$ at larger $(\mu_{q}/T)_{c}$. For nonzero vector repulsion, this reduction is weakened, owing to the shift of the CEP to lower temperature. We note, that for $(\mu_{q}/T)_{c}<0.5$, the $\chi_{B}^{1,2}$ ratio is hardly modified by the repulsive interactions. This can be traced back to a cancellation between the suppression of $\chi_{B}^{2}$ at nonzero $g_{\omega}$ (19) and the shift of the chemical potential (14); both axes are scaled by the factor $(1+9T^{2}G_{\omega}\tilde{\chi}_{B}^{2}(T,\nu))$, for small $\mu_{q}$. On the other hand, the ratio $\chi_{B}^{3}/\chi_{B}^{2}$ is reduced by the factor $(1+9T^{2}G_{\omega}\tilde{\chi}_{B}^{2}(T,\nu))$, for nonvanishing vector interaction. Hence, this ratio is reduced at all values of $\mu_{q}/T$.

As seen in Fig. 4, the skewness ratio, $\chi_{B}^{3,1}$, is an increasing function of $(\mu_{q}/T)_{c}$ along the phase boundary and at the CEP it diverges 222Since the CEP is a singular point, the corresponding value of $\chi_{B}^{3}$ depends on how the limit is taken. Owing to the curvature of the phase boundary, the CEP is approached slightly from below, which implies that $\chi_{B}^{3}$ diverges to plus infinity.. There is also a clear suppression of this ratio along the phase boundary due to repulsive interactions, as seen also for $\chi_{B}^{3,2}$. At small $(\mu_{q}/T)_{c}=0$, the skewness and kurtosis ratios are equal to each other, and differ very little up to $(\mu_{q}/T)_{c}\simeq 0.5$. For larger $(\mu_{q}/T)_{c}$, the kurtosis is decreasing and skewness increasing along the phase boundary, in agreement with recent LGT results Karsch (2017). This behavior also reflects their properties near the critical point. As the CEP is approached along the phase boundary, the kurtosis diverges to minus infinity, while the skewness diverges to plus infinity, as noted above. Thus, both ratios becomes less singular as $g_{\omega}$ is increased, i.e., as the CEP is shifted to larger $\mu_{q}$.

The characteristics of the various net-baryon-number susceptibilities on the phase boundary, shown in Fig. 4, are expected to be similar in QCD. This is because, they are, by and large, determined by $O(4)$ critical fluctuations and by ”confinement”, which are both common to QCD and the PQM model. This opens the possibility to verify directly, if these features of criticality are also reflected in the data on net-proton-number fluctuations obtained in heavy-ion collisions by the STAR Collaboration at several RHIC energies Adamczyk et al. (2014); Luo (2015, 2016).

Clearly, a direct comparison of model results with data has to be taken with caution. Although, the model provides a viable description of the dynamics that drives the system towards chiral symmetry restoration, the spectrum of hadronic degrees of freedom in the low temperature phase is incomplete. Moreover, net-baryon-number fluctuations are in nucleus-nucleus collision experiments quantified by the net proton number. It has been extensively discussed, to what extent are net-proton-number fluctuations accurate proxyies for those of the net baryon number Bzdak and Koch (2012); Bzdak et al. (2013); Kitazawa and Asakawa (2012a, b). Furthermore, there are kinematical cuts, imposed on the STAR data on the cumulants of net proton number, which are not accounted for in the model results. However, these differences are, to a large extent, eliminated by considering ratios of susceptibilities. This assumption is supported by the behavior of the ratio $\chi_{B}^{1,2}$, which, in spite of the differences discussed above, is well approximated by $\tanh(3\mu_{q}/T)$ (see Fig. 5) in LQCD, in the PQM model as well as in the STAR data. The fact that $\chi_{B}^{1,2}$ is well approximated by this functional form indicates that in the transition region, the effective degrees of freedom with nonvanishing baryon number have $B=\pm 1$.

A comparison of results obtained in the PQM model with data requires a correspondence between the collision energy $\sqrt{s}$ and the thermal parameters $(\mu_{q},T)$. Here we employ the phenomenological relation, obtained by analysing the freeze-out conditions in terms of the hadron-resonance-gas model (HRG) Braun-Munzinger et al. ; Andronic et al. (2006, 2009); Cleymans et al. (2006). We then use the resulting dependence of $\mu_{B}$ and $T$ on $\sqrt{s}$ to assign a value for the ratio $(\mu_{q}/T)$ to each of the STAR beam energies. We note that, for $\mu_{q}/T<1$, the phenomenological freeze-out line coincides within errors with the crossover phase boundary obtained in lattice QCD Borsanyi et al. (2010); Andronic et al. (2009). This motivates 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 in Ref. Bazavov et al. (2012b).

In Fig. 5, we show the STAR data on net-proton-number susceptibility ratios and the corresponding PQM model results on net-baryon-number fluctuations computed along the phase boundary. The model results for the ratios $\chi_{B}^{1,2}$, $\chi_{B}^{3,1}$ and $\chi_{B}^{3,2}$ are in qualitative agreement with the data in the whole energy range. For the kurtosis ratio, $\chi_{B}^{4,2}$, this is the case also up to the SPS energy, i.e., for $\sqrt{s}\geq 20$ GeV. However, for $\mu_{q}/T>0.5$, the data on the kurtosis ratio exhibits a qualitatively different dependence on $\mu_{q}/T$ than expected for the critical behavior of $\chi_{B}^{4,2}$, as the CEP is approached along the phase boundary.

As noted above, the ratio $\chi_{B}^{1,2}$ is, on the phase boundary, approximately given by $\chi_{B}^{1,2}\simeq\tanh(3\mu_{q}/T)$ up to $\mu_{q}/T\simeq 1$, as seen in Fig. 5. This form of the ratio of the lowest order cumulants is also obtained in the HRG model and is consistent with LQCD Karsch and Redlich (2011); Allton et al. (2005). Fig. 5 reveals that the data are consistent with this form as well. This fact clearly supports the use of the $\sqrt{s}$ dependence of $\mu_{q}/T$ obtained from the HRG model analysis of multiplicities in the comparison of the PQM model results with data.

In order to assess the sensitivity of the fluctuation observables to model parameters, we show in Fig. 5 ratios of cumulants for four sets of model parameters. The different sets are obtained by varying the sigma meson mass and the form of the Polyakov-loop potential. The parameter sets are described in Appendix B. Fig. 5 shows that, although some quantitative differences can be identified at larger $\mu_{q}/T$, the skewness, kurtosis and $\chi_{B}^{1,2}$ ratios along the phase boundary are qualitatively similar for the different sets of model parameters.

In the comparison of model predictions with data in Fig. 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 understanding of the data. In order to obtain a more quantitative description, we follow Refs. Bazavov et al. (2012b) and Almási et al. (2016); Fu and Pawlowski (2016), and determine the freeze-out conditions by fitting the data on the $\chi_{B}^{3,1}$ ratio, using the $(\sqrt{s})$-dependence of $\mu_{q}/T$ obtained from the fit of the HRG model to particle multiplicities Andronic et al. (2006, 2009); Cleymans et al. (2006).

In Fig. 6 we show the fluctuation ratios along the freeze-out line, which is fixed through the skewness data. The model results are obtained for the four sets of initial conditions. Fig. 6 clearly shows that, along the freeze-out line, the spread of all fluctuations ratios considered for the various parameter sets is much weaker than that observed in Fig. 5 along the phase boundary. This indicates, that moderate changes of the sigma mass and modifications of the form of the Polyakov loop potential may lead to a shift in the temperature scale but essentially with no change of the relative structure of the cumulant ratios.

The results presented in Fig. 6 clearly show that the model provides a very good description of the data on $\chi_{B}^{1,2}$ and $\chi_{B}^{3,2}$. Also the kurtosis data, obtained at higher collision energies, are consistent with model results. However, at $\sqrt{s}<20$ GeV they exhibit a different trend, with the data increasing rapidly at lower energies, while the model result keeps decreasing. We conclude that an increase of $\chi_{B}^{4,2}$ ratio beyond unity, observed in the STAR data at $\sqrt{s}<20$ GeV, is not expected in equilibrium on the chiral critical line nor on the freeze-out line.

As noted above, the ratio $\chi_{B}^{6,2}$ is particularly interesting for identifying criticality governed by the $O(4)$ universality class. This is seen in Fig. 1 at $\mu_{q}=0$, where $\chi_{B}^{6,2}$ is negative at $T=T_{c}$. At $\mu_{q}>0$, the influence of criticality is even more pronounced, as shown in Fig. 7. There, $\chi_{B}^{6,2}$ exhibits a highly nonmonotonic structure near the phase boundary.

In the left panel of Fig. 7 we show the temperature dependence of the ratio $\chi_{B}^{6,2}$ near the phase boundary for $\mu_{q}=0.1$ GeV. This ratio exhibits a very strong variation near the transition temperature. It is deeply negative below $T_{c}$ and develops a maximum just above $T_{c}$. The value of $\chi_{B}^{6,2}$ is thus very sensitive to the freeze-out temperature.

On the right of Fig. 7 we show the ratios $\chi_{B}^{4}/\chi_{B}^{2}$, $\chi_{B}^{5}/\chi_{B}^{1}$ and $\chi_{B}^{6}/\chi_{B}^{2}$ computed, using the model parameters of set A, along the freeze-out line which is determined by fitting $\chi_{B}^{3}/\chi_{B}^{1}$ to data. At vanishing $\mu_{q}$, $\chi_{B}^{4}/\chi_{B}^{2}=\chi_{B}^{3}/\chi_{B}^{1}$ and $\chi_{B}^{6}/\chi_{B}^{2}=\chi_{B}^{5}/\chi_{B}^{1}$, as a consequence of Eq. (22), while at larger $\mu_{q}/T$, these ratios separate.

We note at this point, that the equality of the ratios $\chi_{B}^{4}/\chi_{B}^{2}$ and $\chi_{B}^{3}/\chi_{B}^{1}$ in the STAR data at the highest energy is a strong indication that the fluctuations probed by these cumulants are in thermal equilibrium. It is very unlikely that a system not in equilibrium would yield ratios of cumulants that satisfy Eq. (22). Note, that at $\mu_{q}=0$, the critical $O(4)$ fluctuations yield divergent contributions only to $\chi^{6}_{B}$ and higher cumulants Friman et al. (2011). Thus, the fluctuations probed by $\chi_{B}^{4}/\chi_{B}^{2}$ and $\chi_{B}^{3}/\chi_{B}^{1}$ are not critical. However, a measurement of the ratios involving the fifth and sixth order cumulants would probe whether the $O(4)$ critical fluctuations are in equilibrium or not. Obviously, this test of equilibration is meaningful only at small $\mu_{q}/T$, i.e., only at the highest energies.

At moderate values of $\mu_{q}/T$, the $\chi_{B}^{6}/\chi_{B}^{2}$ ratio is negative and deviates clearly from $\chi_{B}^{5}/\chi_{B}^{1}$. At still lower energies, it exhibits a strong increase towards the CEP, where it diverges. Similarly, $\chi_{B}^{5}/\chi_{B}^{1}$ decreases at moderate $\mu_{q}/T$, and increases strongly as the CEP is approached. These results indicate that in heavy ion collisions $\chi_{B}^{6}/\chi_{B}^{2}$ and $\chi_{B}^{5}/\chi_{B}^{1}$ will exhibit strong nonmonotonic dependencies on $\sqrt{s}$.

Recently, first results on $\chi_{B}^{6}/\chi_{B}^{2}$ in Au-Au collisions at $\sqrt{s}=200\;\mathrm{GeV}$ where reported by the STAR Collaboration for several centralities Esha . The data show a strong suppression of fluctuations compared to the kurtosis ratio. In mid-central and the most central collisions, the $\chi_{B}^{6}/\chi_{B}^{2}$ fluctuation ratio is negative, albeit with still very large statistical uncertainties. In fact, given the large errors for the most central collisions, the preliminary data is consistent with a vanishing $\chi_{B}^{6}$. A value close to zero is consistent with the model results, shown in the right panel of Fig. 7.

The comparison of model results on ratios of net-baryon-number susceptibilities with the STAR data in Figs. 5, 6, and 7, shows that the data, with the exception of kurtosis at low energies, follow general trends expected due to critical chiral dynamics. We note, that the ratios of net-baryon-number susceptibilities near the phase boundary involving net-baryon number cumulants $\chi_{B}^{n}$ with $n\geq 3$ are controlled mainly by the scaling functions in the $O(4)$ and $Z(2)$ universality classes, respectively. This observation indicates, that by measuring fluctuations of conserved charges in heavy-ion collisions, we are indeed probing the QCD phase boundary, and thus accumulating evidence for chiral symmetry restoration.

However, as discussed above, there are several uncertainties and assumptions which must be thoroughly understood before the QCD phase boundary can be pinned down with confidence. Possible contributions to fluctuation observables from effects not related with critical phenomena, like e.g. baryon-number conservation Bzdak et al. (2013) and volume fluctuations Skokov et al. (2013); Braun-Munzinger et al. (2017); Palhares et al. (2011); Fraga et al. (2011); Palhares and Fraga (2012); Hippert et al. (2016) are being explored. We note, however, that volume fluctuations do not modify the double ratios proposed as tests of equilibration at $\mu_{q}\to 0$. We also mention the rather strong sensitivity of higher order net-proton-number cummulants on the transverse momentum range imposed in the analysis of the STAR data. Nevertheless, it is intriguing, that the dynamics of the model, provides a good description of the STAR data (except for $\chi_{B}^{4}$ at the lowest energies), without all these effects of non-critical origin. It remains an important task to assess the effect of theses additional sources of fluctuations in the whole energy range probed by the expriments.

IV Summary and Conclussions

We have studied the influence of the chiral phase transition on fluctuation observables in strongly interacting medium at finite temperature and net baryon number density. We focused on the properties of net-baryon-number fluctuations which are quantified by the $n^{\mathrm{th}}$-order susceptibilities, $\chi_{B}^{n}$. The cumulants $\chi_{B}^{n}$ are directly influenced by the chiral phase transition due to the coupling of the quarks to the scalar sigma field. Furthermore, since the $\chi_{B}^{n}$ are accessible experimentally, they are ideal observables to identify the phase boundary and the critical structures in the QCD phase diagram.

The dynamics was modelled with the Polyakov loop extended Quark-Meson (PQM) model. To correctly account for the critical behavior at the chiral symmetry restoration transition in the $O(4)$ and $Z(2)$ universality classes, we employed the Functional Renormalization Group (FRG). We formulated the FRG equations in the presence of repulsive interactions, and derived the flow equations for derivatives of the thermodynamic pressure.

The main point of our studies was to identify the relations between $\chi_{B}^{n}$ susceptibilities in different $(T,\mu_{B})$ - regions of the phase diagram. To reduce the influence of the non-critical characteristics of the model, like e.g. the mass spectrum or the kinematical cuts on particle momentum distributions, we have computed ratios of susceptibilities. Here, of particular interest, are ratios of the first and second order cumulants, $\chi_{B}^{1}/\chi_{B}^{2}$, the skewness, $\chi_{B}^{3}/\chi_{B}^{1}$, the kurtosis, $\chi_{B}^{4}/\chi_{B}^{2}$, and the sixth-order cumulant $\chi_{B}^{6}/\chi_{B}^{2}$. The higher-order ($n\geq 3$) cumulants are probes of chiral criticality and can be reconstructed experimentally in nucleus-nucleus collisions.

We have calculated, for the first time, the contour plots in the $(T,\mu_{B})$-plane of the above ratios in the PQM model within FRG approach. We have quantified the systematic relations along the phase boundary and the phenomenological freeze-out line, extracted by fitting the skewness data. The influence of the repulsive interactions and the position of the critical endpoint on different $\chi_{B}^{n}$ ratios, was also explored for the first time.

The model results were confronted with ratios of cumulants of net proton number measured by the STAR Collaboration in nucleus-nucleus collisions at energies ranging from $\sqrt{s}=7.7$ GeV up to 200 GeV.

Considering that the relations between different susceptibilities along the phase boundary are induced in part by the universal scaling functions common to QCD and PQM models, and the phenomenological observation that freeze-out in heavy-ion collisions appears near the QCD phase boundary, we have compared the experimental data and model results. The main objective was to verify if the systematics between net-baryon-number susceptibilities observed in the model calculations along the phase boundary are also reflected in the experimental data. To compare model predictions and data we have assumed that the $\sqrt{s}$-dependence of $\mu_{q}/T$, follows the phenomenological freeze-out conditions in nucleu-nucleus collisions. The validity of this assumption is supported by the observation that, for $\mu_{q}/T<1$, the ratio of the first and second order cumulants is well approximated by $\chi_{B}^{1}/\chi_{B}^{2}=\tanh(3\mu_{q}/T)$, in the model, in LQCD and in the data, using the adopted relation between beam energy and $\mu_{q}/T$.

We have shown, that the STAR data for the ratios composed of $\chi_{B}^{n}$ with $n=1,2,3$ follow generic expectations, and thus that they are consistent with criticality at the chiral phase boundary. The quantitative agreement between the calculated cumulant ratios and STAR data is improved when the freeze-out line is determined by fitting $\chi_{B}^{3}/\chi_{B}^{1}$ to data. The data on the kurtosis ratio are also consistent with model systematics at energies ${\sqrt{s}\geq 20\;\mathrm{GeV}}$. However, the strong enhancement of $\kappa\sigma^{2}>1$, observed in nucleus-nucleus collisions at lower energies, is not reproduced by the model along the path in the $(\mu_{q}/T)$- plane, where the ratios of lower order cumulants, $\chi_{B}^{n}$ with $n=1,2,3$, are well described. We have shown, that this conclusion is not affected by the initial conditions nor by the parametrization of the Polyakov-loop potential.

Clearly, the comparison of model results and data is biased by the various assumptions and uncertainties, discussed in Sec. III. However, in spite of all uncertainties, the ratios of net-proton-number susceptibilities obtained by STAR are, as shown in this study, largely consistent with the systematics expected near the chiral phase boundary.

It is interesting to note that a recent analysis of the STAR data on skewness and kurtosis ratios, at $\sqrt{s}\geq 19.6$ GeV, using LQCD results for the Taylor series in $\mu_{q}/T$, shows that the data exhibits all features expected in QCD near the phase boundary Karsch (2017). This result support our findings in PQM model that near chiral phase boundary, ratios of net-baryon-number cumulants, capture properties, which are common to QCD and the PQM model.

Acknowledgments

We acknowledge stimulating discussions with Peter Braun-Munzinger, Frithjof Karsch and Nu Xu. G.A. is also grateful to Vladimir Skokov for discussions on vector interactions. The work of B.F. and K.R. was partly supported by the Extreme Matter Institute EMMI. K. R. also acknowledges partial supports of the Polish National Science Center (NCN) under Maestro grant DEC-2013/10/A/ST2/00106. G. A. acknowledges the support of the Hessian LOEWE initiative through the Helmholtz International 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 $\chi_{B}^{1}$ and $\chi_{B}^{2}$ cumulants of net-baryon-number fluctuations. We start from the flow equation (II.1) for the grand canonical potential,

[TABLE]

where the boson and fermion contribution is introduced, as

[TABLE]

with the pion and sigma masses obtained from

[TABLE]

The flow equations for derivatives, introduced in Eq. (II.1), are derived for each gridpoint from Eq. (23), as

[TABLE]

where $\lambda_{i}$ and $\lambda_{j}$ stand for one of the $\{\mu,\Phi,\bar{\Phi}\}$ parameters, and the prime denotes derivative with respect to $m^{2}$.

Appendix B Polyakov-loop potential and initial conditions

To verify possible model dependence of our results, we have considered different parameterizations of the Polyakov-loop potential.

We apply the polynomial Polyakov-loop potential Ratti et al. (2006), as

[TABLE]

where $b_{3}=0.75$, $b_{4}=7.5$, and

[TABLE]

with $a_{0}=6.75$, $a_{1}=-1.95$, $a_{2}=2.625$, $a_{3}=-7.44$, and $T_{0}=270\,\mathrm{MeV}$.

The parameter $T_{0}$ corresponds to the value of the critical temperature at the first order deconfinement phase transition in a pure Yang-Mills theory. Effects of unquenching can be taken into account by changing $T_{0}$. In the case of two quark flavors, we use $T_{0}=208\,\textrm{MeV}$ from Ref. Herbst et al. (2011).

The second potential applied in our calculations, accounts for the color group structure, and reproduces lattice results on the equation of state and fluctuations of the Polyakov loops, in a pure SU(3) gauge theory. The potential is parameterized Lo et al. (2013), as

[TABLE]

where

[TABLE]

The temperature dependent coefficients are given by

[TABLE]

where $t=T/T_{0}$. The numerical values of the constants in Eq. (32), are summarized in Tab. 1. We set ${T_{0}=270\;\textrm{MeV}}$ for this potential. Lower values of $T_{0}$ would result in distinct position of the peaks of the Polyakov loop and chiral susceptibilities, which is in contrast with LQCD results.

The initial conditions for the flow equations are determined by the vacuum conditions:

[TABLE]

The Yukawa coupling $g_{s}$ and the external field are then determined from

[TABLE]

However, the values of $\lambda$ and $m^{2}$, in the initial meson potentials, are still undefined. They can be fixed by specifying the cutoff scale $\Lambda$, where the flow is started, and the mass of the sigma meson, $m_{\sigma}$. We apply, two distinct values of the initial cutoff scales, and the sigma masses, as

[TABLE]

where the parameters in Eq. (36), were introduced in Ref. Herbst et al. (2013).

From the two parametrization of the Polyakov-loop potentials, Eqs. (28) and (30), and the two chiral initial conditions defined in Eqs. (35) and (36), we define in Table 2, the four different parameter sets that are used in our calculations to study the influence of the model assumptions on the final results. If the parameter set is not explicitly specified, then results refer to set A.