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

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.


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 [1,2] at the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Laboratory and the future Compressed Baryonic Matter (CBM) experiment [3] 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 [4][5][6][7][8][9]. The corresponding pseudocritical temperature has been localized at T c ≈ 156 MeV within a definition-dependent range of several MeV [6,7,10,11]. 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 [12][13][14][15][16][17][18][19][20][21]; see, e.g., Refs. [22,23] for review articles. This notion is supported by results from Dyson-Schwinger equations [24][25][26][27][28], see Ref. [29] 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 freezeout 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 [13,14,[30][31][32][33][34]. 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. [35,36] for reviews. Lattice-QCD results for fluctuations and correlations at zero [37][38][39][40][41] and small chemical potential [42,43] are available, but need to be extended toward higher chemical potential.
A different approach is possible in functional methods. In a series of works [25][26][27]56] a coupled system of Dyson-Schwinger equations (DSEs) for the quark and gluon propagators has been considered and the physics of the Columbia plot [57] has been explored. Results for QCD with heavy quarks and at physical quark masses (N f = 2 + 1 and N f = 2 + 1 + 1) but zero chemical potential agree with corresponding lattice results, see Ref. [29] for an overview. A critical endpoint has been found at T CEP , µ CEP B = (117, 488) MeV that corresponds to a ratio µ CEP B / T CEP ≈ 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 µ B = 0 to and beyond the CEP. We thereby improve previous results for fluctuations calculated in the DSE framework of Ref. [58], 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 [25,56]. In Sec. IV, we present our results and finally conclude in Sec. V.

II. FLUCTUATIONS
In N f = 2 + 1 flavor QCD, there is a conserved charge for each quark flavor controlled by the three quark chemical potentials µ u , µ d , and µ s . The quantities under study in the present work are fluctuations of these conserved charges, i.e., higher-order derivatives of the grandcanonical potential with respect to the quark chemical potentials. Here, Z is the partition function of QCD, T the temperature, and V the volume of the system. The fluctuations are then included in 1 with i, j, k ∈ N 0 . The quark chemical potentials are related to the ones for baryon number (B), strangeness (S), and electric charge (Q) via With these relations one finds for example the secondorder baryon number fluctuation (6) 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 eventby-event analyses of heavy-ion collisions, see the review 1 We 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., χ uds ijk ≡ χ uds ijk (T, µu, µ d , µs). If a subscript is vanishing, it is omitted together with its superscript counterpart, e.g., χ u 2 ≡ χ uds 200 .
articles [35,36,59] for more details. Interesting ratios related to the baryon number are where κ B , σ 2 B , S B , and M 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 [13, 14, 31-34, 60, 61]. At the critical endpoint, the correlation length ξ diverges (at least for infinite volume) and χ B 2 ∼ ξ c with c > 0. From the first BES at RHIC, the STAR collaboration extracted results for the net-proton number fluctuations M P , σ 2 P , S P , and κ P [62,63], 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 [42,43,64] 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 ≤ i + j + k ≤ 4 and determine them via the quark number densities. We start with the grandcanonical potential expressed as a functional of the propagators of QCD. Consequently, Ω contains contributions from quarks, gluons, and ghosts. Since the latter are only weakly chemical-potential dependent, we neglect their contributions and arrive at [65] 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 Φ 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 δ Ω/δS = 0. The quark number densities then read 2 where f ∈ {u, d, s} labels the flavor, N c = 3 denotes the number of colors, Z f 2 is the quark wave function renormalization constant, and q = (ω q , q) with fermionic Matsubara frequencies ω q = (2 q + 1)πT ; q ∈ Z. The Matsubara sum as well as the three-momentum integration is abbreviated by q ≡ T q ∈ Z d 3 q /(2π) 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 = ψ † ψ f = ψ γ 4 ψ f , and its gauge invariance as a density of a global charge is guaranteed by the Landau-Khalatnikov-Fradkin transformations [66,67].
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. [28,68] that is an Euclidean version of the contour-integration technique for Matsubara sums [69]. The regularized quark number density is given by The last term does not depend explicitly on temperature or chemical potential and is known as a "vacuum contribution" in the literature [69]. 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 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 [58,70] 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 [71,72], to N f = 2 [25,73,74], and finally to N f = 2 + 1 and N f = 2 + 1 + 1 quark flavors with physical quark masses [25,26].
With O(4) symmetry broken to O(3) due to the heat bath, the dressed inverse quark propagator S −1 f for a flavor f at nonzero temperature T and quark chemical potential µ f is given by with momentum p = (ω p , p) and fermionic Matsubara frequencies 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 with Z f m denoting the quark mass renormalization constant, and m f is the renormalized current quark mass. In principle, there is a fourth Dirac structure γ 4 γ · p contributing to the inverse quark propagator. However, its contribution is negligible [75] 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 = (ω k , k), where ω k = 2 k πT ( k ∈ 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 with projectors The dressed quark and the dressed gluon propagator each satisfy a Dyson-Schwinger equation that read The symbol D YM µν denotes the sum of the inverse bare gluon propagator and all diagrams with no explicit quark content. The quark self-energy Σ f and the gluon self- energy from the quark loop Π µν are given by with f ∈ {u, d, s}. The equations are shown diagrammatically in Figs. 1 and 2. Note that different flavors are nontrivially coupled through the quark loop Π µν . Furthermore, k = q − p and p = q − k in the quark and gluon DSE, respectively, Z 3 is the ghost renormalization constant, and Γ f ν denotes the dressed quark-gluon vertex. For the coupling we use α = g 2 /(4π) = 0.3 (see Refs. [71,72] for details) and C F = (N 2 c − 1)/(2N c ) is the SU(N 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 YM µν , by quenched temperature-dependent lattice data [72,76]. This approximation misses implicit quark-loop effects in the Yang-Mills self-energies, which are subleading in a 1/N c expansion as compared to the explicit quark loop Π µν . At zero temperature, the effects of this approximation can be estimated using the framework of Ref. [77] and are found to be well below the five-percent level. This strategy has been used in Refs. [26,27,56] 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. [26] for more details): The leading term of the Ball-Chiu vertex construction [78] is multiplied by a phenomenological vertex dressing function Γ that accounts for non-Abelian effects and the correct logarithmic running of the propagators in the ultraviolet. The resulting equations are and The squared momentum argument of Γ is x = k 2 in the quark selfenergy but x = p 2 +q 2 in the quark loop. This is necessary to maintain multiplicative renormalizability of the gluon DSE [77]. 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 GeV 2 and Λ = 1.4 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 u,d,s are fixed using lattice results for the subtracted quark condensate. The quark condensate is given by 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 FIG. 3. Left: Subtracted quark condensate normalized to its vacuum value as a function of temperature at vanishing chemical potential compared to the continuum-extrapolated lattice result of Ref. [6]. Right: Our result for the phase diagram for N f = 2 + 1 quark flavors compared to freeze-out points from heavy-ion collisions extracted by different methods/groups [79][80][81][82][83][84]. Shown is also the region of the chiral crossover from lattice QCD (blue band) [10] (see also Ref. [11]).
up-to-strange quark mass ratio of m s /m u = 25.7 using results for the pion and kaon masses in vacuum obtained from the Bethe-Salpeter formalism developed in Ref. [85]. This results in m u (ζ) = 0.8 MeV and m s (ζ) = 20.56 MeV at a renormalization point of ζ = 80 GeV, i.e., far in the perturbative regime. Note that the values of d 1 and m u,s are slightly different from the ones reported in Ref. [26] because we employ a slightly lighter strange quark mass and an Pauli-Villars regulator 3 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 µ Q and µ S on µ B . For temperatures around 150 MeV, the leading-order result from lattice QCD is µ Q ≈ −0.02 µ B while µ S ≈ 0.2 µ B [40,86]. Thus, to a good approximation we choose µ u = µ d . Furthermore, it has been checked within the framework of DSEs that values for the strange quark chemical potential between µ s = 0 and µ s = µ d hardly affects the location of the CEP [87]. Thus, for the purpose of this work we choose µ s = 0 and, the baryon chemical potential is then given by µ B = 3µ u . 4 3 This amounts to Dµν (k) → Dµν (k)/(1 + k 2 /Λ 2 PV ) in Eq. (20). For the Pauli-Villars scale we use Λ PV = 200 GeV. 4 In Ref. [88], 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.

A. Phase diagram
Before we discuss fluctuations, we present our updated result for the QCD phase diagram with N f = 2 + 1 quark flavors, which closely resembles the one already published in Refs. [26,29]. 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., and find 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 [25] with lattice results [89], as discussed and summarized in Ref. [29]. 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 secondorder CEP at and [38], respectively. Right: Up-quark number density in the vicinity of the CEP for three different chemical potentials.
followed by the coexistence region (shaded gray) of a first-order transition bound by spinodals (solid black). 5 Furthermore, we show the line of baryon chemical potential to temperature ratio µ B /T = 3 (dotted black), emphasizing that the CEP occures at rather large chemical potential with a ratio of µ CEP B / T CEP ≈ 4.2. Our updated value for the location of the critical endpoint is only slightly different than the previous DSE result of Ref. [26]. Again, the error in Eq. (28) 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 [91]. 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) [10] (see also Ref. [11]). 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 µ B /T ≈ 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 µ B /T ≤ 2−2.5 [10,92] in agreement with our result. Furthermore, we also show results for the freezeout points extracted from heavy-ion collisions by different groups/methods [79][80][81][82][83][84]. In the crossover region at small chemical potential, the different results for the freezeout points spread over almost 20 MeV in temperature at and around the pseudocritical temperature extracted on the lattice (see [40,86] 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 firstorder 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 [27].

B. 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 [7,38]. 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 ≈ 165 MeV and reaches an asymptotic value of approximately 0.72 in the high-temperature region. Clearly, this saturation is below the Stefan-Boltzmann limit (χ u 2 → 1 as T → ∞) and happens at much too low temperatures. We attribute this to a known deficiency of our quark-gluon interaction. 6 Most 6 Due to numerical efficiency, our vertex ansatz Eq. (22) takes only the leading Dirac structure γν 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 Tc. 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 important for the purpose of this work, however, is the temperature region 120 -160 MeV below and around the crossover temperature, where the vertex ansatz delivers satisfying results.
Since χ u 2 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 (χ2) c = 153 MeV, which is only slightly lower than the value from the inflection point of the subtracted condensate determined in the previous section to 156 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 [6,7,10,11,95]. While we do not expect χ u 2 to be strongly affected by our choice of µ s = 0 (instead of implementing strangeness neutrality), this is certainly different for χ s 2 . 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 u = T 3 χ u 1 as a function of temperature for three different chemical potentials around our CEP, Eq. (28). For µ u = µ CEP B / 3 = 165 MeV (dashed red), the slope tends to infinity at T = T 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. [96,97]). The location of the phase transition lies within the (at this chemical potential principle, this behavior could be mimicked by making the vertex strength parameter temperature dependent, i.e., d 1 = d 1 (T ) [93]. 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. [87,94] and need to be corroborated. 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 µ u < µ CEP u . In general, at large temperatures the density is an almost linear function of the temperature regardless of the value of the chemical potential.

C. 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 offdiagonal elements that are much harder to be determined and are relegated to future work. 7 Then, the n th -order baryon number fluctuation is given by with n ≥ 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 µ B = 0 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. 8 The behavior of χ B 2 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, χ B 2 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 χ B 3 /χ B 2 (left) and the kurtosis ratio χ B 4 /χ 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 8 While 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.
decrease in magnitude when the temperature is further increased. This structure becomes extremely pronounced close to the CEP. Correspondingly, the kurtosis ratio χ B 4 /χ 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 [35,36] 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 χ B 1 /χ B 2 extracted along our crossover transition line. For small chemical potential, i.e., up to µ B /T 1.5, it is expected from the HRG model [98,99] that the ratio is approximately given by tanh(µ B /T ). This has been seen as well in the PQM model [100] and also shows up in our calculation. Sizeable deviations only occur for larger chemical potential: After the maximum at µ B /T ≈ 1.6, the ratio goes down again and signals the approach to the CEP due the increase of χ B 2 already seen in Fig. 5. Even more interesting are the ratios involving higherorder fluctuations. In Fig. 8 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 √ s = 14.5 GeV on, about halfway toward our CEP, this agreement becomes worse and disappears for √ s ≤ 11.5 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 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 √ s = 19.6 GeV and √ s = 14.5 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 √ s > 14.5 GeV. The results for the kurtosis ratio at √ s = 11.5 GeV and √ s = 7.7 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 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 (bottom) along the crossover line and for lines with a fixed temperature distance from the crossover. Also shown are preliminary data from the STAR collaboration [63,101] at most central collisions. We adopt the µ B -√ s translation from Ref. [82].
already above and are reviewed in Refs. [35,36]. Our theoretical calculation suffers from several limitations. First, we did not yet take into account the effect of offdiagonal 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 [27,94] and/or systematic comparisons with similar calculations in the functional renormalization group framework [91,102,103]. 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. [73], it has been shown that the inclusion of fluctuations from composite pion and sigma fields in the quark DSE serves to generate the critical 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 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 [97], 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. [90].

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 f = 2 + 1 quark flavors, which has been studied extensively in the past [29] and shown to agree with lattice results in the small and moderate chemicalpotential region. At large chemical potential, where lattice QCD cannot be applied, this approach features a critical endpoint at T CEP , µ CEP B = (119, 495) 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 µ B < 250 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.