Baryon number fluctuations in the QCD phase diagram from Dyson-Schwinger equations

TL;DR

This paper investigates baryon number fluctuations in QCD at finite temperature and chemical potential using Dyson-Schwinger equations, providing insights into the critical endpoint and comparing with experimental data.

Contribution

It applies Dyson-Schwinger equations to compute baryon fluctuations in QCD, exploring their behavior near the critical endpoint and comparing with STAR experimental results.

Findings

Fluctuations increase near the critical endpoint.

Results are compatible with a critical endpoint at high chemical potential.

Skewness and kurtosis ratios align with experimental data.

Abstract

We present results for fluctuations of the baryon number for QCD at nonzero temperature and chemical potential. These are extracted from solutions to a coupled set of truncated Dyson-Schwinger equations for the quark and gluon propagators of Landau-gauge QCD with $N_f = 2 + 1$ quark flavors, that has been studied previously. We discuss the changes of fluctuations and ratios thereof up to fourth order for several temperatures and baryon chemical potential up to and beyond the critical endpoint. In the context of preliminary STAR data for the skewness and kurtosis ratios, the results are compatible with the scenario of a critical endpoint at large chemical potential and slightly offset from the freeze-out line. We also discuss the caveats involved in this comparison.

Figures

Equations

$$\Omega =-\frac{T}{V}\log \mathcal{Z}(T,{\mu }_{\text{u}},{\mu }_{\text{d}},{\mu }_{\text{s}})$$

$${\chi }_{ijk}^{\text{u}\text{d}\text{s}}=-\frac{1}{T^{4-(i+j+k)}}\frac{{\partial }^{i+j+k}\Omega }{\partial {\mu }_{\text{u}}^i\partial {\mu }_{\text{d}}^j\partial {\mu }_{\text{s}}^k}$$

$${\mu }_{\text{u}}$$

$${\mu }_{\text{d}}$$

$${\mu }_{\text{s}}$$

$${\chi }_2^{\text{B}}$$

$$=\frac{1}{9}\left[{\chi }_2^{\text{u}}+{\chi }_2^{\text{d}}+{\chi }_2^{\text{s}}+2\left({\chi }_{11}^{\text{u}\text{s}}+{\chi }_{11}^{\text{d}\text{s}}+{\chi }_{11}^{\text{u}\text{d}}\right)\right]$$

$$\begin{array}{c}\frac{{\chi }_4^{\text{B}}}{{\chi }_2^{\text{B}}}={\kappa }_{\text{B}}{\sigma }_{\text{B}}^2,\frac{{\chi }_3^{\text{B}}}{{\chi }_2^{\text{B}}}=S_{\text{B}}{\sigma }_{\text{B}},\\ \frac{{\chi }_1^{\text{B}}}{{\chi }_2^{\text{B}}}=\frac{M_{\text{B}}}{{\sigma }_{\text{B}}^2},\end{array}$$

$$\Omega=-\frac{T}{V}\left(\operatorname{Tr}\log\frac{S^{-1}}{T}-\operatorname{Tr}\bigl{[}\openone-S_{0}^{-1}S\bigr{]}+\Phi_{\textup{int}}[S]\right)$$

$$n_{f}=-\frac{\partial\hskip 0.83344pt\Omega}{\partial\mu_{f}}=-N_{\textup{c}}\hskip 0.83344ptZ_{2}^{f}\hskip 0.83344pt\operatorname{\mathchoice{\ooalign{\raisebox{0.08pt}{\scalebox{0.85}{$\displaystyle\sum$}}\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.24pt}{\scalebox{0.6}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.24pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.24pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{q}\operatorname{Tr}\bigl{[}\gamma_{4}\hskip 0.83344ptS_{f}(q)\bigr{]}\,,$$

$$n_f^{\text{reg}}=-N_{\text{c}}{Z}_2^f\int \frac{{\text{d}}^3\mathbf{q}}{(2\pi {)}^3}{K}_f(\mathbf{q})$$

$$K_f(\mathbf{q})$$

$$\displaystyle\phantom{=\;}-\frac{1}{2\pi}\int_{-\infty}^{\infty}\textup{d}q_{4}\operatorname{Tr}\bigl{[}\gamma_{4}\hskip 0.83344ptS_{f}(q_{4},\bm{q})\bigr{]}\,.$$

$${\chi }_2^{\text{u}}=\frac{1}{T^2}\frac{\partial n_{\text{u}}^{\text{reg}}}{\partial {\mu }_{\text{u}}}$$

$$S_f^{-1}(p)=\text{i}({\omega }_p+\text{i}{\mu }_f){\gamma }_4{C}_f(p)+\text{i}\mathbf{\gamma }\cdot \mathbf{p}{A}_f(p)+B_f(p)$$

$$S_{0,f}^{-1}(p)=Z_{2}^{f}\hskip 0.83344pt\bigl{(}\hskip 0.83344pt\textup{i}\hskip 0.83344pt(\omega_{p}+\textup{i}\hskip 0.83344pt\mu_{f})\hskip 0.83344pt\gamma_{4}+\textup{i}\hskip 0.83344pt\bm{\gamma}\cdot\bm{p}\hskip 0.83344pt+Z_{m}^{f}\hskip 0.83344ptm_{f}\hskip 0.83344pt\bigr{)}$$

$$D_{\mu \nu }(k)=P_{\mu \nu }^{\text{T}}(k)\frac{Z_{\text{T}}(k)}{k^2}+P_{\mu \nu }^{\text{L}}(k)\frac{Z_{\text{L}}(k)}{k^2},$$

$$P_{\mu \nu }^{\text{T}}(k)$$

$$P_{\mu \nu }^{\text{L}}(k)$$

$$S_f^{-1}(p)=S_{0,f}^{-1}(p)+{\Sigma }_f(p),$$

$$\displaystyle D_{\mu\nu}^{-1}(k)=\bigl{[}D_{\mu\nu}^{\textup{YM}}(k)\bigr{]}^{-1}+\Pi_{\mu\nu}(k)\,.$$

$${\Sigma }_f(p)$$

$${\Pi }_{\mu \nu }(k)$$

$$\frac{1}{{\tilde{Z}}_3}{\Gamma }_{\nu }^f(q,p;k)$$

$$\displaystyle\phantom{=\;}+(1-\delta_{4\nu})\hskip 0.83344pt\frac{A_{f}(q)+A_{f}(p)}{2}\hskip 0.83344pt\biggr{]}\hskip 0.83344pt\gamma_{\nu}$$

$$\Gamma (x)=\frac{d_1}{d_2+x}+\frac{1}{1+{\Lambda }^2/x}{\left(\frac{\alpha {\beta }_0}{4\pi }\log \left(1+\frac{x}{{\Lambda }^2}\right)\right)}^{2\delta },$$

$$\langle\bar{\psi}\psi\rangle_{f}=-N_{\textup{c}}\hskip 0.83344ptZ_{2}^{f}Z_{m}^{f}\hskip 0.83344pt\operatorname{\mathchoice{\ooalign{\raisebox{0.08pt}{\scalebox{0.85}{$\displaystyle\sum$}}\cr$\displaystyle\int$\cr}}{\ooalign{\raisebox{0.24pt}{\scalebox{0.6}{$\textstyle\sum$}}\cr$\textstyle\int$\cr}}{\ooalign{\raisebox{0.24pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}{\ooalign{\raisebox{0.24pt}{\scalebox{0.6}{$\scriptstyle\sum$}}\cr$\scriptstyle\int$\cr}}}_{q}\operatorname{Tr}\bigl{[}S_{f}(q)\bigr{]}$$

$${\Delta }_{\text{u}\text{s}}=⟨\bar{\psi }\psi {⟩}_{\text{u}}-\frac{m_{\text{u}}}{m_{\text{s}}}⟨\bar{\psi }\psi {⟩}_{\text{s}},$$

$$T_{\text{c}}=\underset{T}{\mathrm{argmax}}\left|\frac{\partial {\Delta }_{\text{u}\text{s}}}{\partial T}\right|$$

$$T_{\textup{c}}=$156\pm 1.\text{\,}\mathrm{MeV}$$

$

$$\bigl{(}T^{\textup{CEP}},\,\mu_{\textup{B}}^{\textup{CEP}}\bigr{)}=($119\pm 2.\text{\,}$,$495\pm 2.\text{\,}$)\,$\mathrm{MeV}$$

$

$$\chi_{n}^{\textup{B}}=\frac{1}{3^{n}}\,\bigl{(}2\hskip 0.83344pt\chi_{n}^{\textup{u}}+\chi_{n}^{\textup{s}}\bigr{)}$$

Full text

I Introduction

Extracting the location of a putative critical endpoint (CEP) of QCD from heavy-ion collisions is one of the major goals of the Beam Energy Scan (BES) program Aggarwal et al. (2010a); STAR collaboration (2014) at the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory and the future Compressed Baryonic Matter (CBM) experiment Friman et al. (2011a) at the Facility for Antiproton and Ion Research (FAIR).

Theoretically it is by no means clear that such a critical endpoint exists. At zero chemical potential, there is firm evidence from lattice QCD for an analytic crossover from a low-temperature phase characterized by chiral symmetry breaking to a high-temperature (partially) chirally restored phase Aoki et al. (2006, 2009); Borsányi et al. (2010); Bazavov et al. (2012a); Bhattacharya et al. (2014); Bazavov et al. (2014). The corresponding pseudocritical temperature has been localized at $T_{\textup{c}}\approx$156\text{,}\mathrm{MeV}$$ within a definition-dependent range of several MeV Borsányi et al. (2010); Bazavov et al. (2012a); Bellwied et al. (2015a); Bazavov et al. (2019). However, the situation is much less clear at (real) chemical potential, where lattice calculations are hampered by the notorious fermion sign problem. Model calculations suggest that the continuous crossover becomes steeper with increasing chemical potential and finally merges into a second-order CEP followed by a region of a first-order phase transition at large chemical potential Asakawa and Yazaki (1989); Stephanov et al. (1998, 1999); Fukushima (2004); Megias et al. (2006); Ratti et al. (2006); Schaefer et al. (2007); Skokov et al. (2010a, b); Herbst et al. (2011); see, e.g., Refs. Drews and Weise (2017); Fukushima and Skokov (2017) for review articles. This notion is supported by results from Dyson-Schwinger equations Qin et al. (2011); Fischer and Luecker (2013); Fischer et al. (2014); Eichmann et al. (2016); Gao and Liu (2016), see Ref. Fischer (2019) for a recent review.

In order to put these theoretical ideas to the test in experiments, observables have been identified that are capable to deliver signals of the CEP. Provided the freeze-out in heavy-ion collisions is sufficiently close to the phase boundary, fluctuations of conserved charges (baryon number, strangeness, and electric charge) are expected to provide this information Stephanov et al. (1998, 1999); Jeon and Koch (1999, 2000); Asakawa et al. (2000); Koch et al. (2005); Ejiri et al. (2006). Various ratios of cumulants of these conserved quantities can be extracted from experiment in event-by-event analyses and compared to corresponding ratios of fluctuations that can be determined in theoretical calculations, see, e.g., Refs. Luo and Xu (2017); Bzdak et al. (2019) for reviews. Lattice-QCD results for fluctuations and correlations at zero Cheng et al. (2009); Borsányi et al. (2012); Bazavov et al. (2012b); Borsányi et al. (2013); Bellwied et al. (2015b) and small chemical potential Bazavov et al. (2017a); Borsányi et al. (2018) are available, but need to be extended toward higher chemical potential.

In hadron resonance gas (HRG) approaches and refined effective models, such as the Polyakov-loop enhanced Nambu-Jona-Lasinio model and the Polyakov-loop quark-meson model (PQM), a wealth of interesting results on fluctuations have been obtained already, see, e.g., Refs. Skokov et al. (2010b, 2011); Karsch et al. (2011); Schaefer and Wagner (2012); Morita et al. (2015); Morita and Redlich (2015); Fu and Pawlowski (2015, 2016); Fu et al. (2016); Shao et al. (2018); Fu et al. (2018a); Szymański et al. (2019); Yin et al. (2019) and references therein. Concerning the Yang-Mills sector, these models rely on the Polyakov loop potential that couples aspects of confinement to the chiral dynamics, however without backcoupling. Thus, gluons are no active degrees of freedom (d.o.f.) and their reaction to the medium can neither be studied nor directly taken into account.

A different approach is possible in functional methods. In a series of works Fischer and Luecker (2013); Fischer et al. (2014, 2015); Eichmann et al. (2016) a coupled system of Dyson–Schwinger equations (DSEs) for the quark and gluon propagators has been considered and the physics of the Columbia plot Brown et al. (1990) has been explored. Results for QCD with heavy quarks and at physical quark masses ($N_{\textup{f}}=2+1$ and $N_{\textup{f}}=2+1+1$) but zero chemical potential agree with corresponding lattice results, see Ref. Fischer (2019) for an overview. A critical endpoint has been found at \bigl{(}T^{\textup{CEP}},\,\mu_{\textup{B}}^{\textup{CEP}}\bigr{)}=(117,488)\,\mathrm{MeV}$$ that corresponds to a ratio $\mu_{\textup{B}}^{\textup{CEP}}/\,T^{\textup{CEP}}\approx 4.2$, i.e., large chemical potential. In this work we will use this framework to explore cumulants and ratios thereof along the crossover line from $\mu_{\textup{B}}=0$ to and beyond the CEP. We thereby improve previous results for fluctuations calculated in the DSE framework of Ref. Xin et al. (2014), where backcoupling effects have not been taken into account. In particular, we discuss ratios involving the skewness and kurtosis and compare our results with preliminary data from the STAR collaboration extracted from the BES at RHIC.

This work is organized as follows: In Sec. II, we detail our method to extract fluctuations from the quark propagator and derivatives thereof. In Sec. III, we summarize the truncation scheme of the DSEs and discuss the (slight) changes as compared to previous works Fischer and Luecker (2013); Fischer et al. (2015). In Sec. IV, we present our results and finally conclude in Sec. V.

II fluctuations

In $N_{\textup{f}}=2+1$ flavor QCD, there is a conserved charge for each quark flavor controlled by the three quark chemical potentials $\mu_{\textup{u}}$, $\mu_{\textup{d}}$, and $\mu_{\textup{s}}$. The quantities under study in the present work are fluctuations of these conserved charges, i.e., higher-order derivatives of the grand-canonical potential

[TABLE]

with respect to the quark chemical potentials. Here, $\mathcal{Z}$ is the partition function of QCD, $T$ the temperature, and $V$ the volume of the system. The fluctuations are then included in111We usually suppress the arguments of the fluctuations. However, one has to keep in mind that they are functions of temperature and all chemical potentials, i.e., $\chi_{ijk}^{\textup{u}\textup{d}\textup{s}}\equiv\chi_{ijk}^{\textup{u}\textup{d}\textup{s}}(T,\mu_{\textup{u}},\mu_{\textup{d}},\mu_{\textup{s}})$. If a subscript is vanishing, it is omitted together with its superscript counterpart, e.g., $\chi_{2}^{\textup{u}}\equiv\chi_{200}^{\textup{u}\textup{d}\textup{s}}$.

[TABLE]

with $i,\hskip 0.83344ptj,\hskip 0.83344ptk\in\mathbb{N}_{0}$. The quark chemical potentials are related to the ones for baryon number (B), strangeness (S), and electric charge (Q) via

[TABLE]

With these relations one finds for example the second-order baryon number fluctuation

[TABLE]

in terms of quark degrees of freedom. Other fluctuations can be determined analogously.

Ratios of fluctuations in baryon number, electric charge, and strangeness are particularly interesting since they are equal to corresponding ratios of cumulants that can be extracted from experimental quantities accessible in event-by-event analyses of heavy-ion collisions, see the review articles Asakawa and Kitazawa (2016); Luo and Xu (2017); Bzdak et al. (2019) for more details. Interesting ratios related to the baryon number are

[TABLE]

where $\kappa_{\textup{B}}$, $\sigma_{\textup{B}}^{2}$, $S_{\textup{B}}$, and $M_{\textup{B}}$ denote the kurtosis, variance, skewness, and mean of the net-baryon distribution, respectively. Ratios of fluctuations are a suitable tool to explore the phase diagram of QCD since they are sensitive to phase transitions Stephanov et al. (1998, 1999); Asakawa et al. (2000); Jeon and Koch (2000); Stephanov (2004); Koch et al. (2005); Ejiri et al. (2006); Friman et al. (2011b). At the critical endpoint, the correlation length $\xi$ diverges (at least for infinite volume) and $\chi_{2}^{\textup{B}}\sim\xi^{c}$ with $c>0$.

From the first BES at RHIC, the STAR collaboration extracted results for the net-proton number fluctuations $M_{\textup{P}}$, $\sigma^{2}_{\textup{P}}$, $S_{\textup{P}}$, and $\kappa_{\textup{P}}$ Aggarwal et al. (2010b); Luo (2015), which can be used as a proxy for fluctuations of the net-baryon number. The data suggest a number of interesting tendencies that are drastically different from results of HRG model calculations, but agree with results from lattice QCD D’Elia et al. (2017); Bazavov et al. (2017a); Borsányi et al. (2018) obtained for small baryon chemical potential. As already mentioned in the introduction, it is the purpose of this paper to provide theoretical results for larger chemical potential in the framework of functional approaches to QCD.

In the present work, we consider the lowest-order fluctuations with $1\leq i+j+k\leq 4$ and determine them via the quark number densities. We start with the grand-canonical potential expressed as a functional of the propagators of QCD. Consequently, $\Omega$ contains contributions from quarks, gluons, and ghosts. Since the latter are only weakly chemical-potential dependent, we neglect their contributions and arrive at Cornwall et al. (1974)

[TABLE]

with $S$ being the dressed quark propagator with flavor, color, Dirac, and momentum degrees of freedom; $S_{0}$ denotes its bare counterpart. The trace has to be taken in the functional sense over flavor, color, Dirac, and momentum space, and the interaction functional $\Phi_{\textup{int}}$ contains all two-particle irreducible diagrams with respect to $S$. Note that Eq. (8) is nothing but the two-particle irreducible effective action evaluated at the stationary point and therefore $\delta\hskip 0.83344pt\Omega\hskip 0.83344pt/\hskip 0.83344pt\delta S=0$.

The quark number densities then read222We work in four-dimensional Euclidean space-time with Hermitian gamma matrices obeying $\{\gamma_{\mu},\gamma_{\nu}\}=2\hskip 0.83344pt\delta_{\mu\nu}$; $\mu,\nu\in\{1,2,3,4\}$. Our choice for the heat bath vector is $u=(u_{4},\bm{u})=(1,\bm{0})$.

[TABLE]

where $f\in\{\textup{u},\textup{d},\textup{s}\}$ labels the flavor, $N_{\textup{c}}=3$ denotes the number of colors, $Z_{2}^{f}$ is the quark wave function renormalization constant, and $q=(\omega_{q},\bm{q})$ with fermionic Matsubara frequencies $\omega_{q}=(2\hskip 0.83344pt\ell_{q}+1)\hskip 0.83344pt\pi T$; $\ell_{q}\in\mathbb{Z}$. The Matsubara sum as well as the three-momentum integration is abbreviated by \operatorname{\mathchoice{\ooalign{\raisebox{0.08pt}{\scalebox{0.85}{\displaystyle\sum}}\cr\displaystyle\int\cr}}{\ooalign{\raisebox{0.24pt}{\scalebox{0.6}{\textstyle\sum}}\cr\textstyle\int\cr}}{\ooalign{\raisebox{0.24pt}{\scalebox{0.6}{\scriptstyle\sum}}\cr\scriptstyle\int\cr}}{\ooalign{\raisebox{0.24pt}{\scalebox{0.6}{\scriptstyle\sum}}\cr\scriptstyle\int\cr}}}_{q}\equiv T\sum_{\ell_{q}\in\hskip 0.5834pt\mathbb{Z}}\int\textup{d}^{3}\bm{q}\hskip 0.83344pt/\hskip 0.83344pt(2\pi)^{3}, and the remaining trace has to be evaluated in Dirac space. The quantity $S_{f}(q)$ denotes the dressed quark propagator at nonzero temperature and chemical potential with only Dirac and momentum degrees of freedom left. Equation (9) can also be written as an expectation value, $n_{f}=\langle\psi^{\dagger}\psi\rangle_{f}=\langle\bar{\psi}\hskip 0.83344pt\gamma_{4}\hskip 0.83344pt\psi\rangle_{f}$, and its gauge invariance as a density of a global charge is guaranteed by the Landau–Khalatnikov–Fradkin transformations Landau and Khalatnikov (1956); *Fradkin:1955jr; Zumino (1960).

The quark number densities need to be evaluated using a very large number of Matsubara frequencies in order to obtain stable results. In addition, with nonperturbative propagators and a numerical cutoff in the three-momentum integral, this expression needs to be regularized. To this end, we employ the subtraction scheme used in Refs. Gao et al. (2016); Gao and Liu (2016) that is an Euclidean version of the contour-integration technique for Matsubara sums Kapusta and Gale (2006). The regularized quark number density is given by

[TABLE]

with

[TABLE]

The last term does not depend explicitly on temperature or chemical potential and is known as a “vacuum contribution” in the literature Kapusta and Gale (2006). We verified that its subtraction leads to cutoff-independent results.

Having the quark number densities at hand, the fluctuations are obtained by higher-order derivatives of these densities. For example, the second-order up quark fluctuation reads

[TABLE]

and other quantities are obtained analogously.

III Dyson–Schwinger equations

In the following, we briefly summarize the functional framework used to determine the temperature and chemical potential dependent dressed quark propagator needed to determine the quark densities of Eq. (9). To this end, we solve a set of truncated Dyson–Schwinger equations. In contrast to previous works on fluctuations in the DSE framework Xin et al. (2014); Xu et al. (2015) we take the back reaction of the quarks onto the Yang–Mills sector explicitly into account. This establishes a temperature and chemical-potential dependence of the gluon controlled by QCD dynamics rather than simple modelling. Furthermore, this allows for explicit control over the quark-flavor dependence of all results. Our framework evolved from the quenched case Fischer (2009); Fischer et al. (2010), to $N_{\textup{f}}=2$ Fischer and Mueller (2011); Fischer et al. (2011); Fischer and Luecker (2013), and finally to $N_{\textup{f}}=2+1$ and $N_{\textup{f}}=2+1+1$ quark flavors with physical quark masses Fischer and Luecker (2013); Fischer et al. (2014).

With $\textup{O}(4)$ symmetry broken to $\textup{O}(3)$ due to the heat bath, the dressed inverse quark propagator $S_{f}^{-1}$ for a flavor $f$ at nonzero temperature $T$ and quark chemical potential $\mu_{f}$ is given by

[TABLE]

with momentum $p=(\omega_{p},\bm{p})$ and fermionic Matsubara frequencies $\omega_{p}=(2\hskip 0.83344pt\ell_{p}+1)\hskip 0.83344pt\pi T$, $\ell_{p}\in\mathbb{Z}$. The dressing functions $C_{f}$, $A_{f}$, and $B_{f}$ depend on momentum and contain all nonperturbative information. Furthermore, they depend on temperature and chemical potential. The corresponding bare quark propagator reads

[TABLE]

with $Z_{m}^{f}$ denoting the quark mass renormalization constant, and $m_{f}$ is the renormalized current quark mass. In principle, there is a fourth Dirac structure $\gamma_{4}\hskip 0.83344pt\bm{\gamma}\cdot\bm{p}$ contributing to the inverse quark propagator. However, its contribution is negligible Contant and Huber (2017) and therefore not considered in this work.

Since we work in Landau gauge, the gluon propagator is purely transverse with respect to its four-momentum $k=(\omega_{k},\bm{k})$, where $\omega_{k}=2\hskip 0.83344pt\ell_{k}\pi T$ ($\ell_{k}\in\mathbb{Z}$) are bosonic Matsubara frequencies. However, due to the presence of the heat bath, the transverse space splits into two parts, and the dressed gluon propagator reads

[TABLE]

with projectors

[TABLE]

The dressed quark and the dressed gluon propagator each satisfy a Dyson–Schwinger equation that read

[TABLE]

The symbol $D_{\mu\nu}^{\textup{YM}}$ denotes the sum of the inverse bare gluon propagator and all diagrams with no explicit quark content. The quark self-energy $\Sigma_{f}$ and the gluon self-energy from the quark loop $\Pi_{\mu\nu}$ are given by

[TABLE]

with $f\in\{\textup{u},\textup{d},\textup{s}\}$. The equations are shown diagrammatically in Figs. 2 and 2. Note that different flavors are nontrivially coupled through the quark loop $\Pi_{\mu\nu}$. Furthermore, $k=q-p$ and $p=q-k$ in the quark and gluon DSE, respectively, $\widetilde{Z}_{3}$ is the ghost renormalization constant, and $\Gamma_{\nu}^{f}$ denotes the dressed quark-gluon vertex. For the coupling we use $\alpha=g^{2}/\hskip 0.83344pt(4\pi)=0.3$ (see Refs. Fischer (2009); Fischer et al. (2010) for details) and $C_{\textup{F}}=(N_{\textup{c}}^{2}-1)\hskip 0.83344pt/\hskip 0.83344pt(2N_{\textup{c}})$ is the $\textup{SU}(N_{\textup{c}})$ quadratic Casimir factor in the fundamental representation that stems from the color trace.

In order to determine the gluon propagator at nonzero temperature and chemical potential, an efficient approximation that has been used in the literature is to replace the Yang–Mills part of the equation, $D_{\mu\nu}^{\textup{YM}}$, by quenched temperature-dependent lattice data Fischer et al. (2010); Maas et al. (2012). This approximation misses implicit quark-loop effects in the Yang–Mills self-energies, which are subleading in a $1\hskip 0.83344pt/\hskip 0.83344ptN_{\textup{c}}$ expansion as compared to the explicit quark loop $\Pi_{\mu\nu}$. At zero temperature, the effects of this approximation can be estimated using the framework of Ref. Fischer and Alkofer (2003) and are found to be well below the five-percent level. This strategy has been used in Refs. Fischer et al. (2014, 2015); Eichmann et al. (2016) to determine the location of the critical endpoint and will be adopted also in this work.

Eventually, the quark-gluon vertex is the last quantity that needs to be specified to obtain a closed system of equations. We use the following ansatz (see Ref. Fischer et al. (2014) for more details): The leading term of the Ball–Chiu vertex construction Ball and Chiu (1980) is multiplied by a phenomenological vertex dressing function $\Gamma$ that accounts for non-Abelian effects and the correct logarithmic running of the propagators in the ultraviolet. The resulting equations are

[TABLE]

and

[TABLE]

where $\delta=-9N_{\textup{c}}\hskip 0.83344pt/\hskip 0.83344pt(44N_{\textup{c}}-8N_{\textup{f}})$ is the anomalous dimension of the vertex and $\beta_{0}=(11N_{\textup{c}}-2N_{\textup{f}})\hskip 0.83344pt/\hskip 0.83344pt3$. The squared momentum argument of $\Gamma$ is $x=k^{2}$ in the quark self-energy but $x=p^{2}+q^{2}$ in the quark loop. This is necessary to maintain multiplicative renormalizability of the gluon DSE Fischer and Alkofer (2003). Note that the vertex ansatz includes effects from nonzero temperature and chemical potential, as the full vertex certainly would, since the terms from the Ball–Chiu construction involve the quark dressing functions. The parameters $d_{2}=$0.5\text{,}{\mathrm{GeV}}^{2} and $\Lambda=$1.4\text{\,}\mathrm{GeV} are fixed to match the scales in the quenched gluon propagator from the lattice.

Quark masses, vertex strength, and chemical potentials

The remaining value of the vertex strength parameter $d_{1}$ as well as the quark masses $m_{\textup{u},\textup{d},\textup{s}}$ are fixed using lattice results for the subtracted quark condensate. The quark condensate is given by

[TABLE]

for each quark flavor $f$. It is plagued by a quadratic divergence for all flavors with nonzero quark masses and needs to be regularized. This can be accomplished by the difference

[TABLE]

which defines the subtracted quark condensate. It is an order parameter for chiral symmetry breaking and may be used to define the pseudocritical temperature, see Eq. (42) below. We adapt the vertex strength parameter $d_{1}=$8.49\text{,}{\mathrm{GeV}}^{2} such that the pseudocritical temperature found on the lattice Borsányi *et al.* ([2010](#bib.bib6)); Bazavov *et al.* ([2012a](#bib.bib7)); Bellwied *et al.* ([2015a](#bib.bib10)); Bazavov *et al.* ([2019](#bib.bib11)) is reproduced. We work in the isospin-symmetric limit of equal up and down quark masses, $m_{\textup{u}}=m_{\textup{d}}$. In the high-temperature phase, $\Delta_{\textup{u}\textup{s}}$ is mainly controlled by these masses and we adapt their values to match the lattice results. Finally we fix the up-to-strange quark mass ratio of $m_{\textup{s}}\hskip 0.83344pt/\hskip 0.83344ptm_{\textup{u}}=25.7$ using results for the pion and kaon masses in vacuum obtained from the Bethe–Salpeter formalism developed in Ref. Heupel *et al.* ([2014](#bib.bib86)). This results in $m_{\textup{u}}(\zeta)=$0.8\text{\,}\mathrm{MeV} and $m_{\textup{s}}(\zeta)=$20.56\text{,}\mathrm{MeV} at a renormalization point of $\zeta=$80\text{\,}\mathrm{GeV}, i.e., far in the perturbative regime. Note that the values of $d_{1}$ and $m_{\textup{u},\textup{s}}$ are slightly different from the ones reported in Ref. Fischer et al. (2014) because we employ a slightly lighter strange quark mass and an Pauli–Villars regulator333This amounts to $D_{\mu\nu}(k)\to D_{\mu\nu}(k)\hskip 0.83344pt/\hskip 0.83344pt(1+k^{2}/\hskip 0.83344pt\Lambda_{\textup{PV}}^{2})$ in Eq. (28). For the Pauli–Villars scale we use $\Lambda_{\textup{PV}}=$200\text{,}\mathrm{GeV}$$.

for the quark DSE rather than a hard cutoff.

In principle, the chemical potentials should be adjusted in order to implement strangeness neutrality as encountered in a heavy-ion collision. This is done by an appropriate dependence of $\mu_{\textup{Q}}$ and $\mu_{\textup{S}}$ on $\mu_{\textup{B}}$. For temperatures around $150\text{\,}\mathrm{MeV}$, the leading-order result from lattice QCD is $\mu_{\textup{Q}}\approx-0.02\,\mu_{\textup{B}}$ while $\mu_{\textup{S}}\approx 0.2\,\mu_{\textup{B}}$ Borsányi et al. (2013); Bazavov et al. (2012c). Thus, to a good approximation we choose $\mu_{\textup{u}}=\mu_{\textup{d}}$. Furthermore, it has been checked within the framework of DSEs that values for the strange quark chemical potential between $\mu_{\textup{s}}=0$ and $\mu_{\textup{s}}=\mu_{\textup{d}}$ hardly affects the location of the CEP Welzbacher (2016). Thus, for the purpose of this work we choose $\mu_{\textup{s}}=0$ and, the baryon chemical potential is then given by $\mu_{\textup{B}}=3\hskip 0.83344pt\mu_{\textup{u}}$.444In Ref. Fu et al. (2018b), the impact of strangeness neutrality on thermodynamic observables is studied. Since at large chemical potentials quantitative corrections of the order of 20% for some thermodynamic quantities have been found, which may also affect fluctuations, we strive to implement strangeness neutrality in future work.

IV Results

IV.1 Phase diagram

Before we discuss fluctuations, we present our updated result for the QCD phase diagram with $N_{\textup{f}}=2+1$ quark flavors, which closely resembles the one already published in Refs. Fischer et al. (2014); Fischer (2019). Overall, changes due to the slightly adapted strange quark mass and the different regularization scheme are very small.

We determine the pseudocritical temperature of the chiral crossover from the inflection point of the subtracted quark condensate with temperature, i.e.,

[TABLE]

and find

[TABLE]

at vanishing chemical potential. The error given is purely numerical in nature. In the left diagram of Fig. 3, we show the subtracted quark condensate as a function of temperature at vanishing chemical potential. As described in the previous section, our result for the pseudocritical temperature agrees by construction with the lattice result. A nontrivial result, however, is the almost perfect match regarding the steepness of the chiral transition. Another highly nontrivial result is the matching of the unquenched gluon propagator Fischer and Luecker (2013) with lattice results Aouane et al. (2013), as discussed and summarized in Ref. Fischer (2019).

Our result for the phase diagram at nonzero chemical potential is shown in the right diagram of Fig. 3. The chiral crossover line (dashed black) becomes steeper with increasing chemical potential and terminates in a second-order CEP at

[TABLE]

followed by the coexistence region (shaded gray) of a first-order transition bound by spinodals (solid black).555See Ref. Gunkel et al. (2019) for a more detailed discussion of the coexistence region between the spinodal lines. Furthermore, we show the line of baryon chemical potential to temperature ratio $\mu_{\textup{B}}/\hskip 0.83344ptT=3$ (dotted black), emphasizing that the CEP occures at rather large chemical potential with a ratio of $\mu_{\textup{B}}^{\textup{CEP}}/\,T^{\textup{CEP}}\approx 4.2$. Our updated value for the location of the critical endpoint is only slightly different than the previous DSE result of Ref. Fischer et al. (2014). Again, the error in Eq. (44) is purely numerical. In order to estimate the systematic error due to our truncation assumptions, we need to compare with different truncations as, e.g., employed very recently in the framework of functional renormalization-group equations Fu et al. (2019). This will be a task for future work.

In the plot in Fig. 3, we also show results for the chiral transition obtained on the lattice (blue band) Bellwied et al. (2015a) (see also Ref. Bazavov et al. (2019)). As can be seen in the plot, this band features a (very) small error at small chemical potential which rapidly increases toward larger chemical potential. At about $\mu_{\textup{B}}/\hskip 0.83344ptT\approx 3$, the errors become so large that further extrapolation becomes meaningless. Combined evidence of different methods on the lattice points toward no critical endpoint for $\mu_{\textup{B}}/\hskip 0.83344ptT\leq 2-2.5$ Bellwied et al. (2015a); Bazavov et al. (2017b) in agreement with our result. Furthermore, we also show results for the freeze-out points extracted from heavy-ion collisions by different groups/methods Alba et al. (2014); Becattini et al. (2017); Vovchenko et al. (2016); Adamczyk et al. (2017); Andronic et al. (2017, 2018). In the crossover region at small chemical potential, the different results for the freeze-out points spread over almost $20\text{\,}\mathrm{MeV}$ in temperature at and around the pseudocritical temperature extracted on the lattice (see Bazavov et al. (2012c); Borsányi et al. (2013) for a direct comparison of lattice results and experimental data). While in the presence of a (first-order) phase transition one would expect the freeze-out to occur at temperatures below the critical one, this notion is hard to formulate in the crossover region, where no unique definition of a critical temperature exists. At large chemical potential, however, where we see a first-order transition in our DSE results, we have to expect corrections either to the location of the experimental freeze-out points or to the DSE results in order to account for a proper ordering of temperatures. In this respect we would like to point out that potential corrections to the DSE calculation have already been identified (on a qualitative basis), which have the potential to shift the CEP to larger temperatures and/or chemical potentials thereby resolving this problem Eichmann et al. (2016).

IV.2 Quark number fluctuations

After the discussion of the phase diagram in the last subsection, we now focus on our results for the fluctuations. In the left diagram of Fig. 4, the second-order up/down quark fluctuation is shown as a function of temperature at vanishing chemical potential (solid black) and compared to results from lattice QCD Bazavov et al. (2012a); Borsányi et al. (2012). The agreement between both approaches is not as good as for the quark condensate but still very reasonable. The DSE result increases up to temperatures of $T\approx$165\text{,}\mathrm{MeV}$$ and reaches an asymptotic value of approximately $0.72$ in the high-temperature region. Clearly, this saturation is below the Stefan–Boltzmann limit ($\chi_{2}^{\textup{u}}\to 1$ as $T\to\infty$) and happens at much too low temperatures. We attribute this to a known deficiency of our quark-gluon interaction.666Due to numerical efficiency, our vertex ansatz Eq. (38) takes only the leading Dirac structure $\gamma_{\nu}$ into account. In the Landau gauge employed in this work, the full vertex contains 24 different structures. Half of these are only present when chiral symmetry is broken, i.e., these terms react strongly on the chiral restoration around $T_{\textup{c}}$. This effect is not captured by the ansatz. As a result, the continuous weakening of the quark-gluon interaction that drives the system and its fluctuations toward the Stefan–Boltzmann limit is not properly represented, thus leading to the high-temperature artifacts seen in the left diagram of Fig. 4. In principle, this behavior could be mimicked by making the vertex strength parameter temperature dependent, i.e., $d_{1}=d_{1}(T)$ Mueller et al. (2010). Such a modification could be motivated and guided, e.g., by an explicit calculation of (parts of) the vertex as a function of temperature. Preliminary results of this endeavor have been discussed in Refs. Welzbacher (2016); Contant et al. (2018) and need to be corroborated. Most important for the purpose of this work, however, is the temperature region $120$ – $160\text{\,}\mathrm{MeV}$ below and around the crossover temperature, where the vertex ansatz delivers satisfying results.

Since $\chi_{2}^{\textup{u}}$ experiences the most rapid growth in the region of the chiral crossover, its inflection point with temperature can also be used to define the pseudocritical temperature. We find $T_{\textup{c}}^{\,(\chi_{2})}=$153\text{,}\mathrm{MeV}$$, which is only slightly lower than the value from the inflection point of the subtracted condensate determined in the previous section to $156\text{\,}\mathrm{MeV}$. Note again that there is no unique definition of the critical temperature due to the crossover nature of the transition and both values are in agreement with lattice QCD Borsányi et al. (2010); Bazavov et al. (2012a); Endrődi et al. (2011); Bellwied et al. (2015a); Bazavov et al. (2019). While we do not expect $\chi_{2}^{\textup{u}}$ to be strongly affected by our choice of $\mu_{\textup{s}}=0$ (instead of implementing strangeness neutrality), this is certainly different for $\chi_{2}^{\textup{s}}$. We therefore postpone a comparison of this quantity with corresponding lattice results to future work.

Next, we turn to nonzero chemical potential. The right diagram of Fig. 4 displays the up/down quark number density $n_{\textup{u}}=T^{3}\hskip 0.83344pt\chi_{1}^{\textup{u}}$ as a function of temperature for three different chemical potentials around our CEP, Eq. (44). For $\mu_{\textup{u}}=\mu_{\textup{B}}^{\textup{CEP}}/\,3=$165\text{,}\mathrm{MeV}$$ (dashed red), the slope tends to infinity at $T=T^{\textup{CEP}}$ corresponding to a diverging second-order fluctuation. For chemical potentials above the critical one, the system undergoes a first-order phase transition. Thus, the density is discontinuous across the phase boundary and shows a finite jump (solid blue). This behavior is consistent with results obtained in effective models (see, e.g., Refs. Buballa (2005); Schaefer and Wambach (2007)). The location of the phase transition lies within the (at this chemical potential small) region between the upper and lower spinodal line shown in our phase diagram (cf. Fig. 3). Below the critical chemical potential, where the transition is an analytic crossover (dash-dotted green), the slope is finite around the pseudocritical temperature, and the density changes continuously as a function of temperature for all $\mu_{\textup{u}}<\mu_{\textup{u}}^{\textup{CEP}}$. In general, at large temperatures the density is an almost linear function of the temperature regardless of the value of the chemical potential.

IV.3 Baryon number fluctuations

Having the quark number fluctuations at hand, we are now able to compute the baryon number fluctuations. In particular we are interested in the changes induced by growing chemical potential in various ratios of baryon number fluctuations as we approach the critical endpoint. In the present work, we restrict ourselves to quark number fluctuations diagonal in quark flavor and neglect off-diagonal elements that are much harder to be determined and are relegated to future work.777This is justified by lattice results indicating that off-diagonal correlations are subleading as compared to diagonal ones Bellwied et al. (2015b). Then, the $n$th-order baryon number fluctuation is given by

[TABLE]

with $n\geq 1$. We choose a selection of fixed baryon chemical potentials and evaluate the fluctuations as a function of temperature. The results are shown in the left diagram of Fig. 5, where we display the second-order baryon number fluctuation approaching the CEP. At $\mu_{\textup{B}}=$0\text{,}\mathrm{MeV}$$ (solid black), the behavior is similar to the up quark number fluctuation, cf. Fig. 4, i.e., we find a monotonous increase below and around the crossover transition. At nonzero chemical potential about halfway toward the critical endpoint (dashed blue), a bulge begins to develop at and around the pseudocritical transition temperature. This bulge becomes larger as we further increase the chemical potential (dash-dotted red). Close to the location of the CEP, the bulge grows considerably and becomes a sharp peak (dash-dot-dotted green) which finally diverges at the CEP (dotted purple), as expected from the behavior of the quark number density, discussed above.888While in principle one could fine-tune the chemical potential to come arbitrarily close to the actually divergence, in practice limited numerical accuracy together with finite computer resources always lead to a very large but still finite correlation length.

The behavior of $\chi_{2}^{\textup{B}}$ in the first-order region of the phase diagram is shown in the right diagram of Fig. 5. The second-order baryon number fluctuation shows two branches corresponding to the chirally broken solution (solid blue) and partially chirally restored solution (dashed red) of the DSE for the quark propagator. The overlap of the two solutions defines the coexistence region of the first-order transition that is bounded by the spinodals at temperatures indicated by vertical dotted, gray lines. For temperatures above and below the coexistence region, $\chi_{2}^{\textup{B}}$ is only very slowly varying with temperature.

Next, we discuss ratios of fluctuations that are directly related to experimental quantities in heavy-ion collisions through event-by-event analyses (see Eq. (7)). In Fig. 6, we plot the skewness ratio $\chi^{\textup{B}}_{3}/\chi^{\textup{B}}_{2}$ (left) and the kurtosis ratio $\chi^{\textup{B}}_{4}/\chi^{\textup{B}}_{2}$ (right) again as a function of temperature for various lines of constant chemical potential up to the critical endpoint. These show distinctive features. Whereas for small chemical potential up to halfway toward the critical endpoint all structures are very small in size, these grow rapidly when the CEP is approached. The skewness develops a characteristic rise with temperature accompanied by a zero crossing and subsequent equally drastic decrease in magnitude when the temperature is further increased. This structure becomes extremely pronounced close to the CEP. Correspondingly, the kurtosis ratio $\chi^{\textup{B}}_{4}/\chi^{\textup{B}}_{2}$ develops an asymmetric double-peak structure across the phase boundary.

There are a number of caveats when comparing results from theoretical calculations with data extracted from experiment. These are related to the experimental conditions such as the finite volume and the finite temporal extent of the fireball and the question whether and when the system is in thermodynamical equilibrium. Furthermore, these are related to details of the experimental analysis such as centrality cuts, the question whether proton number fluctuations are a proxy for baryon number fluctuations, and potential other issues, see the reviews Luo and Xu (2017); Bzdak et al. (2019) and references therein. Still, there is considerable interest in comparing experimental data with results from theoretical calculations along the phase boundary. Such a comparison is done in the following.

In Fig. 7, we display our results for the ratio $\chi^{\textup{B}}_{1}/\chi^{\textup{B}}_{2}$ extracted along our crossover transition line. For small chemical potential, i.e., up to $\mu_{\textup{B}}/\hskip 0.83344ptT\lesssim 1.5$, it is expected from the HRG model Karsch et al. (2003); Karsch and Redlich (2011) that the ratio is approximately given by $\tanh(\mu_{\textup{B}}/\hskip 0.83344ptT)$. This has been seen as well in the PQM model Almasi et al. (2017) and also shows up in our calculation. Sizeable deviations only occur for larger chemical potential: After the maximum at $\mu_{\textup{B}}/\hskip 0.83344ptT\approx 1.6$, the ratio goes down again and signals the approach to the CEP due the increase of $\chi^{\textup{B}}_{2}$ already seen in Fig. 5.

Even more interesting are the ratios involving higher-order fluctuations. In Fig. 8, we present results for the skewness ratio $\chi^{\textup{B}}_{3}/\chi^{\textup{B}}_{2}$ (upper diagram; blue, solid line) and the kurtosis ratio $\chi^{\textup{B}}_{4}/\chi^{\textup{B}}_{2}$ (lower diagram; blue, solid line) along our chiral phase boundary determined from the inflection point of the chiral condensate, Eq. (42). For the skewness, this criterion leaves us on the left and positive branch of the oscillations shown in Fig. 5. For the kurtosis, however, we probe the (small) negative region around the phase boundary once the chemical potential becomes large. At small chemical potential, there is very good agreement between our results and the (preliminary) data from the STAR collaboration. From $\sqrt{s}=$14.5\text{,}\mathrm{GeV} on, about halfway toward our CEP, this agreement becomes worse and disappears for $\sqrt{s}\leq$11.5\text{\,}\mathrm{GeV}. In order to discuss this aspect further, we also evaluated the skewness and kurtosis ratios on lines with a fixed temperature distance of $3$, $6$, and $9\text{\,}\mathrm{MeV}$ below the crossover line. The general idea of this comparison is to study the impact of two different effects: (i) as mentioned already several times, there is no unique definition of the critical temperature in the crossover region and it is therefore by now means clear, whether a given definition should coincide with the experimental freeze-out line or not; (ii) as the chemical potential becomes larger and a potential CEP is approached, it is also not clear whether the freeze-out line and the crossover line have the same curvature. In other words, it may very well be, that the freeze-out line bends stronger than the crossover line and the distance between the two lines grows with chemical potential.

Taken at face value, our results shown in Fig. 8 seem to support this notion at least on a qualitative level. At small chemical potential, the variations in both ratios with temperature are very small and cannot be discriminated by the data. The two data points at $\sqrt{s}=$19.6\text{,}\mathrm{GeV} and $\sqrt{s}=$14.5\text{\,}\mathrm{GeV}, however, favor a scenario with a freeze-out line very close to the crossover line, and we conclude that this is generally the case for $\sqrt{s}>$14.5\text{,}\mathrm{GeV}. The results for the kurtosis ratio at $\sqrt{s}=$11.5\text{\,}\mathrm{GeV} and $\sqrt{s}=$7.7\text{,}\mathrm{GeV}$$, however, suggest that the freeze-out line in this region of the phase diagram is separated from the crossover line by at least $9\text{\,}\mathrm{MeV}$. The corresponding results for the skewness ratio show the same general trend, although on a less quantitative level than the ones for the kurtosis.

There are several caveats involved in the comparison of the experimental STAR data and our results in Fig. 8. Some caveats on the experimental side have been discussed already above and are reviewed in Refs. Luo and Xu (2017); Bzdak et al. (2019). Our theoretical calculation suffers from several limitations. First, we did not yet take into account the effect of off-diagonal contributions to the baryon number fluctuations. Second, there may be a substantial error associated with the precise location of the critical endpoint. The source of this error is entirely located in the truncation for the quark-gluon vertex and may be reduced in the future by extended DSE calculations Eichmann et al. (2016); Contant et al. (2018) and/or systematic comparisons with similar calculations in the functional renormalization group framework Braun et al. (2011, 2016); Fu et al. (2019). Third, one has to bear in mind that the fluctuations triggering the CEP in this work are gluonic in nature. Consequently, the critical exponents of our CEP are mean field. In Ref. Fischer and Mueller (2011), it has been shown that the inclusion of fluctuations from composite pion and sigma fields in the quark DSE serves to generate the critical $\textup{O}(4)$ physics of the chiral two-flavor theory. We therefore expect that the extension of that framework toward chemical potential places our CEP in the correct $\textup{Z}(2)$ universality class due to the fluctuating sigma field. Furthermore, one may expect a decrease of the size of the critical region around the CEP Schaefer and Wambach (2007), which in turn will drive the results shown in Fig. 8 further toward the STAR data even if the location of the CEP remains unchanged. Since the inclusion of pions and the sigma leads to a significant increase in complexity and CPU time in our calculations, this is left for future work. A first step toward this can be found in Ref. Gunkel et al. (2019).

V Summary and conclusions

In this work, we extracted ratios of cumulants involving the skewness and the kurtosis from baryon number fluctuations at nonzero temperature and chemical potential. To this end, we employed a framework of Dyson–Schwinger equations for $N_{\textup{f}}=2+1$ quark flavors, which has been studied extensively in the past Fischer (2019) and shown to agree with lattice results in the small and moderate chemical-potential region. At large chemical potential, where lattice QCD cannot be applied, this approach features a critical endpoint at \bigl{(}T^{\textup{CEP}},\,\mu_{\textup{B}}^{\textup{CEP}}\bigr{)}=(119,495)\,\mathrm{MeV}. Due to inherent limitations of the truncation scheme used, this value may have systematic errors of at least twenty percent. Future heavy-ion collision experiments such as FAIR/CBM, NICA, and the STAR Fixed-Target program will be able to probe the corresponding region of the QCD phase diagram. In order to facilitate these experiments and to make contact with already existing preliminary data from the BES at RHIC, we determined skewness and kurtosis ratios along our crossover line up to the CEP. Furthermore, we scanned lines of equal temperature distance below the transition line. For chemical potentials $\mu_{\textup{B}}<$250\text{\,}\mathrm{MeV}, our results are in agreement with the STAR data. For larger values we obtain qualitative and quantitative differences when we approach the CEP on the crossover transition line. However, qualitative agreement between our results and the STAR data can be obtained if we assume that the freeze-out line and the transition line separate at larger chemical potential. We also discussed several caveats in this interpretation, which need to be checked in future work.

Acknowledgements.

We thank Bernd-Jochen Schaefer for many stimulating discussions on the physics of fluctuations. Furthermore, we thank Thorsten Steinert, Richard Williams, and Fei Gao for fruitful discussions. This work has been supported by the Helmholtz Graduate School for Hadron and Ion Research for FAIR, the GSI Helmholtzzentrum für Schwerionenforschung, the Helmholtz International Center for FAIR within the LOEWE program of the State of Hesse, and the BMBF under contract 05P18RGFCA. M.B. acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the CRC-TR 211 “Strong-interaction matter under extreme conditions”—project number 315477589—TRR 211. Feynman diagrams were drawn with JaxoDraw Binosi et al. (2009).