Skewness, kurtosis, and the fifth and sixth order cumulants of net baryon-number distributions from lattice QCD confront high-statistics STAR data

We present new results on up to sixth-order cumulants of net baryon-number fluctuations at small values of the baryon chemical potential, μ B , obtained in lattice QCD calculations with physical values of light and strange quark masses. Representing the Taylor expansions of higher-order cumulants in terms of the ratio of the two lowest-order cumulants, M B = σ 2 B ¼ χ B 1 ð T; μ B Þ = χ B 2 ð T; μ B Þ , allows for a parameter-free comparison with data on net proton-number cumulants obtained by the STAR Collaboration in the Beam Energy Scan at RHIC. We show that recent high-statistics data on skewness and kurtosis ratios of net proton-number distributions, obtained at a beam energy ﬃﬃﬃﬃﬃﬃﬃﬃ s NN p ¼ 54 . 4 GeV, agree well with lattice QCD results on cumulants of net baryon-number fluctuations close to the pseudocritical temperature, T pc ð μ B Þ , for the chiral transition in QCD. We also present first results from a next-to-leading-order expansion of fifth- and sixth-order cumulants on the line of the pseudocritical temperatures.


I. INTRODUCTION
The phase diagram of strong-interaction matter at nonzero temperature and nonzero baryon-number density is being explored intensively through numerical calculations performed in the framework of lattice-regularized quantum chromodynamics (QCD) [1], as well as through ultrarelativistic heavy-ion collisions with varying beam energies [2]. At vanishing and small values of the chemical potentials for conserved charges [baryon number (μ B ), electric charge (μ Q ), strangeness (μ S )] it is well established that the transition from the low-temperature hadronic region to the quark-gluon plasma at high temperature is a smooth transition [3] characterized by a pseudocritical temperature, T pc ðμ B Þ [4][5][6][7]. At larger values of the baryon chemical potential it, however, is generally expected that a first-order phase transition line exists, which ends in a second-order critical point [8,9]. This elusive critical point is searched for in the Beam Energy Scan (BES) performed at the Relativistic Heavy Ion collider (RHIC) at Brookhaven National Laboratory [10]. However, its existence as a fundamental property of the theory of strong interactions (QCD) still awaits confirmation.
The pseudocritical line, T pc ðμ B Þ, which distinguishes the low-and high-temperature regimes of strong-interaction matter as described by QCD, has been determined quite accurately in lattice QCD calculations for baryon chemical potentials up to about twice the pseudocritical temperature, μ B ≲ 2T pc ð0Þ ≃ 300 MeV [4][5][6][7]. In our recent analysis we found [7] T pc ðμ B Þ ¼ T 0 with T 0 pc ¼ ð156.5 AE 1.5Þ MeV and κ B 2 ¼ 0.012ð4Þ with a Oðμ 4 B Þ correction that vanishes within errors. At μ B ¼ 0 the pseudocritical temperature turns out to be in good agreement with the freeze-out temperature determined by the ALICE Collaboration at the LHC [11] and the pseudocritical line, T pc ðμ B Þ, is also consistent with freeze-out temperatures determined by the STAR Collaboration during the first BES at RHIC (BES-I) [12], albeit these temperatures have larger statistical errors.
The experimental determination of the freeze-out parameter is based on a measurement of particle yields, i.e., first moments of particle distributions, which in turn are closely related to first-order cumulants of net charge fluctuations. The proximity of freeze-out temperatures and the pseudocritical temperature determined in QCD suggests that the higher-order moments of net charge fluctuations also reflect properties of a thermal medium close to the pseudocritical line. This, however, is not at all well established and many caveats have been discussed suggesting that the relation of higher-order cumulants, measured experimentally, to cumulants of conserved charge fluctuations, calculated in equilibrium QCD thermodynamics, is not at all straightforward [10,13].
Higher-order cumulants of net conserved charge fluctuations are obtained as derivatives of the logarithm of the QCD partition functions with respect to the chemical potentials of conserved charges, ⃗ μ ¼ ðμ B ; μ Q ; μ S Þ, χ X n ðT; ⃗ μÞ ¼ 1 VT 3 ∂ n ln ZðT; ⃗ μÞ ∂μ n X ; X ¼ B; Q; S; ð2Þ withμ ≡ μ=T. These higher-order derivatives become increasingly sensitive to long-range correlations and large fluctuations in the vicinity of a critical point. At least from the theoretical point of view higher-order cumulants thus are ideally suited to search for a possible critical point in the QCD phase diagram [14][15][16]. The BES at RHIC aims at finding evidence for such a critical point through the analysis of e.g., higher-order cumulants of net protonnumber fluctuations which are considered to be good proxies for cumulants of net baryon-number fluctuations. Results, obtained with BES-I at RHIC, indicate that qualitative changes in the behavior of net proton-number fluctuations occur at beam energies ffiffiffiffiffiffiffi ffi s NN p ∼ 20 GeV [17,18]. This may hint at the existence of a critical point for large values of the baryon chemical potential.
While the finding of nonmonotonic behavior of higherorder cumulants of net proton-number fluctuations generated well-justified excitement [17,18], we still need to establish that this behavior is caused by thermal fluctuations in the vicinity of a critical point and that these higher-order cumulants indeed probe thermal conditions at the time of freeze-out. As pointed out in Ref. [19] at least for small values of the baryon chemical potential the first four cumulants of net baryon-number fluctuations, i.e., mean 2 ] are predicted in QCD equilibrium thermodynamics to be related.
This relation, which is only slightly violated in strangeness neutral systems, has been established in lattice QCD calculations using next-to-leading-order (NLO) Taylor expansions of the first four cumulants of net baryon-number fluctuations [19]. The data on cumulants of net-proton number fluctuations, obtained by STAR during BES-I [18] 20] are encouraging. As will be discussed in Sec. IV, these data fulfill the above inequality and the difference of the cumulant ratios given in Eq. (3) agrees with lattice QCD results even on a quantitative level.
We will present here new results on the density dependence of up to sixth-order cumulants of net baryon-number fluctuations. We calculate Taylor series at nonzero values of the baryon-number, electric-charge and strangeness chemical potentials that involve up to eighth-order cumulants. We perform these expansions for the case of strangeness neutral systems, n S ¼ 0, with a ratio of electric charge to baryon number, n Q =n B ¼ 0.4, that is representative of the conditions met in heavy-ion collisions. This allows to construct Taylor expansions for nth-order cumulants, 1 χ B n ðT; μ B Þ, up to Oðμ 8−n B Þ. For the case of the skewness and kurtosis ratios, S B σ 3 B =M B and κ B σ 2 B , respectively, we thus can extend earlier NLO calculations and perform next-to-next-toleading-order (NNLO) expansions that allow to better control truncation effects in the Taylor series. We also present, for the first time, results from NLO calculations for the hyperskewness and hyper-kurtosis (fifth-and sixth-order cumulants) ratios χ B 5 ðT; μ B Þ=χ B 1 ðT; μ B Þ and χ B 6 ðT; μ B Þ=χ B 2 ðT; μ B Þ. We show that these ratios are expected to be negative at ffiffiffiffiffiffiffi ffi s NN p ¼ 54.4 GeV, in contrast to the preliminary findings for sixth-order cumulants of net proton-number fluctuations reported by the STAR Collaboration [20]. 1 Rather than specifying in the argument of χ B n all three chemical potentials, ⃗ μ, we give in the strangeness neutral case only the baryon chemical potential. This paper is organized as follows. In the next section we briefly present our calculational setup, the new statistics collected on lattices of size 32 3 × 8 and 48 3 × 12 and the general fitting ansatz used for fits at fixed values of N τ ¼ 8 and 12, joint fits of these data as well as continuum limit estimates. In Sec. III we present results for Taylor expansions of cumulants of net baryon-number fluctuations that use up to eighth-order cumulants. We compare these results with experimental data for cumulants of net proton-number fluctuations in Sec. IV. Section V contains our conclusions. Explicit expressions for the first four Taylor expansion coefficients of net baryon-number cumulants are given in the Appendix.

II. CALCULATIONAL SETUP
Up to fourth-order cumulants of net baryon-number fluctuations have been calculated previously [19,21,22] in a next-to-leading-order Taylor expansion. In particular, we performed calculations [19] with the highly improved staggered quark [23] discretization scheme for (2 þ 1)flavor QCD with a physical strange quark mass and two degenerate, physical light quark masses. Here we extend these calculations by increasing the number of gauge field configurations generated on lattices of size 32 3 × 8 and 48 3 × 12 by a factor of 3-5 in the transition region and at least a factor of 2 at other values of the temperature. This allows us to calculate up to eighth-order cumulants of net baryon-number, net strangeness and net electric-charge fluctuations, including also their correlations, at vanishing values of the chemical potentials. These cumulants provide expansion coefficients in Taylor series for net baryonnumber cumulants χ B n ðT; ⃗ μÞ. We calculate NLO expansions for fifth-and sixth-order cumulants and obtain NNLO results for third-and fourth-order cumulants. In the case of first-and second-order cumulants, i.e., the mean and variance of net baryon-number distributions, we even obtain next-to-NNLO (NNNLO) results. The set of gauge field ensembles, which has been used in this analysis, and the number of gauge field configurations per ensemble on lattices with temporal extent N τ ¼ 8 and 12 are summarized in Table I.
Results for up to eighth-order diagonal net baryonnumber susceptibilities, χ B n ≡ χ B n ðT; 0Þ, are given in Fig. 1. For the quadratic fluctuations, χ B 2 , we also show results for lattices with temporal extent N τ ¼ 6, which were already used in Ref. [7]. For the eighth-order cumulant, χ B 8 , we only show our results for N τ ¼ 8 as statistical errors on the N τ ¼ 12 data are still too large. The bands shown in these figures give a continuum extrapolation for χ B 2 ðTÞ using data from calculations for three different lattice spacings (aT ¼ 1=N τ ) and a continuum estimate for χ B 4 ðTÞ based on N τ ¼ 8 and 12 data sets. For χ B 6 ðTÞ and χ B 8 ðTÞ we only show spline interpolations of the data obtained on the 32 3 × 8 lattices. Results for these cumulants, obtained from calculations within a noninteracting hadron resonance gas (HRG) model that use resonances from the Particle Data Tables [25] (PDG-HRG) as well as additional resonances calculated within the quark model [26,27] (QM-HRG) are given by lines. The latter list contains additional resonances not (yet) observed experimentally.
We determine the expansion coefficients,χ B;k n ðTÞ, for Taylor series of nth-order cumulants, for the case of vanishing net strangeness density, n S ¼ 0, and an electric-charge to baryon-number ratio, n Q =n B ¼ 0.4. Explicit expressions for the NLO expansion coefficients of up to sixth-order net baryon-number cumulants are given in Ref. [19]. The explicit form of the higherorder expansion coefficients are given in the Appendix. Using the Taylor series for nth-order cumulants, Eq. (4), we construct cumulant ratios with polynomials of order ½k max ; l max , ð5Þ In order to control systematic effects arising from the truncation of the Taylor series expansion for the cumulant ratios R B nm , we calculate these ratios using different orders of the Taylor expansion for the cumulants appearing in the numerator and denominator of these ratios. We analyzed the polynomial ratios for different ½k max ; l max as well as Taylor expansions of the ratios themselves. We find that the TABLE I. Number of gauge field configurations on lattices of size 32 3 × 8 and 48 3 × 12 used in the analysis of up to eighthorder Taylor expansion coefficients. The values of the gauge coupling as well as the strange and light quark mass parameter at these temperature values are taken from Ref. [24], where details on the statistics available on the 24 3 × 6 lattices were also given. All configurations are separated by 10 time units in rational hybrid Monte Carlo simulations [24].
No. of conf. former are more stable at large μ B =T. In the following we will use the ratios of polynomials with ½k max ; l max corresponding to identical orders (LO, NLO, NNLO, NNNLO) of expansions in the cumulants appearing in the numerator and denominator, respectively.
We fit cumulant ratios using a rational polynomial ansatz, where T 0 is some arbitrary scale. When using this rational polynomial ansatz for fits at nonzero μ B we allow for a quadratic μ B dependence of all expansion coefficients, a n ðμ B Þ ¼ a n;0 þ a n;2μ 2 B and similarly for b n ðμ B Þ. When performing joint fits of data on lattices with different sizes and lattice spacings, a, we allow for Oða 2 Þ cutoff corrections that are parametrized in terms of the temporal lattice extent N τ ¼ 1=aT, e.g., with gðT;μ B Þ and hðT;μ B Þ being rational polynomials of the type given in Eq. (6).

A. Mean and variance of net baryon-number fluctuations
We have calculated the ratio of the mean, for systems with vanishing net strangeness, n S ¼ 0, and a net electric-charge to net baryon-number density n Q =n B ¼ 0.4 on lattices with temporal extent N τ ¼ 8 and 12. Using up to eighth-order Taylor expansion coefficients, we can construct Taylor series up to order Oðμ 7 B Þ and Oðμ 6 B Þ for χ B 1 ðT; μ B Þ and χ B 2 ðT; μ B Þ, respectively. Truncating these series at k max and l max ¼ k max − 1, respectively, we construct the ½k max ; l max polynomial ratios which provide leading-order ([1, 0]), next-to-leading-order ( [3,2]) etc., approximations for the ratio of the mean and variance of the distribution for net baryon-number fluctuations, Results for different ½k max ; l max are shown in Fig. 2. The figure shows results obtained on lattices with temporal extent N τ ¼ 8 and 12 at a temperature 2 T ≃ 157 MeV which is close to the pseudocritical temperature at μ B ¼ 0.
We find that cutoff effects are negligible for μ B =T ≤ 1 and remain comparable to the statistical errors for the For further details see text. 2 As is evident from Table I the temperatures differ slightly for the two lattice sizes: T ¼ 156.76 MeV for N τ ¼ 8 and T ¼ 157.13 MeV for N τ ¼ 12, respectively. N τ ¼ 12 data at least up to μ B =T ≃ 1.2. This holds true in the entire temperature range T ∈ ½135 MeV∶175 MeV analyzed by us. Differences in R B 12 constructed from NNLO and NNNLO Taylor series of the cumulants are about 2% for μ B =T ¼ 1.
As the temperature dependence of R B 12 is weak in the temperature range considered by us and also deviations of the μ B dependence from the leading order, linear behavior are moderate we find that using [2,3] rational polynomials in both terms of the fit ansatz given in Eq. (7) are sufficient for obtaining good fits to the data. We performed fits separately for the NNLO and NNNLO data sets at fixed values of T and μ B =T ≤ 1.2. The resulting continuum estimates for R B 12 , evaluated for several values of the temperature in the vicinity of the pseudocritical temperature, T pc ð0Þ, are shown in Fig. 3. We note that the variation with temperature is small. As will be discussed in Sec. IV the results obtained for R B 12 at μ B ≲ 125 MeV are in good agreement with HRG model calculations. For larger values of μ B we find, however, R B;QCD 12 > R B;HRG 12 , which reflects the large deviations of higher-order cumulants, evaluated in QCD at μ B ¼ 0, from the corresponding HRG values.

B. Skewness and kurtosis of net baryon-number fluctuations
While the low-order cumulants Þ and their ratio are in good agreement with HRG model calculations that use noninteracting, point-like hadrons at and below T pc (see also discussion in Sec. IV), this clearly is not the case for higher-order cumulants. This is apparent in calculations of the skewness and kurtosis ratios, which both are unity in noninteracting HRG model calculations, but are known to be significantly smaller in lattice QCD calculations already in the vicinity of the pseudocritical temperature, T pc ð0Þ, at vanishing values of the baryon chemical potential. Moreover, in contrast to the cumulant ratio R B 12 , the ratios R B 31 and R B 42 show a much stronger temperature dependence and a milder dependence on μ B . It thus has been suggested that the ratio R B 12 is well suited to determine the baryon chemical potential from experimental data, while the ratios R B 31 and R B 42 constrain the temperature [21,28].
Using our results for up to eighth-order cumulants of conserved charge fluctuations and correlations, we can construct NNLO expansions for the third-and fourth-order cumulants χ B 3 ðT; μ B Þ and χ B 4 ðT; μ B Þ, where again the electric-charge and strangeness chemical potentials have been fixed by demanding n S ¼ 0 and n Q =n B ¼ 0.4. With this we determine up to NNLO results for the skewness and kurtosis cumulant ratios R B 31 and R B 42 . We again first use our high-statistics data obtained on the N τ ¼ 8 lattices to analyze the effect of truncations of the Taylor expansions at finite orders of μ B . We used the fit ansatz given in Eq. (6)   cover the entire parameter range of relevance for the calculation of these cumulant ratios on the pseudocritical line for μ B =T ≲ 1.
In Ref. [19] we showed that the skewness and kurtosis ratios R B 31 and R B 42 are almost identical at leading order, Oðμ 0 B Þ. The NLO correction to the kurtosis ratio R B 42 , however, is about a factor of 3 larger than that for the skewness ratio R B 31 . Figure 4 suggests that these relations are still well respected by the NNLO results. The slope of R B 42 ðT; μ B Þ as a function ofμ B at fixed T is significantly larger than that of R B 31 ðT; μ B Þ and, in fact, it is still consistent with being about a factor of 3 larger. This is shown in Fig. 5 where we compare the μ B -dependent parts of R B 31 and R B 42 =3. Also shown in this figure are the second derivatives of R B 31 ðT; μ B Þ and R B 42 ðT; μ B Þ=3 with respect to μ B =T.
Compared to the lower-order ratio R B 12 higher-order corrections in the Taylor expansion of R B 31 are significantly larger. In the temperature range shown in Fig. 4  are about a factor of 3 larger.
In Fig. 6 we show results for the skewness and kurtosis ratios R B 31 ðT; μ B Þ and R B 42 ðT; μ B Þ obtained at μ B ¼ 0 on lattices with temporal extent N τ ¼ 8 and 12. Obviously results for N τ ¼ 12 are systematically below those for N τ ¼ 8. This is in accordance with the observed shift of the pseudocritical temperatures [7] to smaller values with increasing N τ or, equivalently, decreasing lattice spacing at fixed temperature aT ¼ 1=N τ . When performing joint fits to the N τ ¼ 8 and 12 data, using the ansatz given in  Eq. (7), we find that within our current statistical errors on the N τ ¼ 12 data we cannot resolve any T or μ B =T dependence of cutoff effects. It thus suffices to use a constant ansatz for the cutoff corrections, i.e., we use gðT; μ B Þ ¼ a 0;0 and a [3,4] rational polynomial for the continuum limit result fðT; μ B Þ. A joint fit to the N τ ¼ 8 and 12 data yields a 0;0 ¼ 3.2ð1.5Þ for R B 31 ðT; μ B Þ and a 0;0 ¼ 3.2ð3.0Þ for R B 42 ðT; μ B Þ. The resulting continuum limit estimates at μ B ¼ 0 are also shown in Fig. 6.
The inset in Fig. 6 (bottom) shows the continuum estimate for the difference R B 42 − R B 31 at μ B ¼ 0 as a function of T. At temperatures below T ≃ 150 MeV this difference is consistent with being zero. In the crossover region, T pc ð0Þ ¼ 156.5ð1.5Þ MeV we find that the difference is slightly positive, R B 42 ðT pc Þ − R B 31 ðT pc Þ ¼ 0.008ð3Þ. Continuum estimates for R B 31 ðT; μ B Þ and R B 42 ðT; μ B Þ at two values of the temperature, corresponding to the current error band for the pseudocritical temperature at μ B ¼ 0 are shown in Fig. 7.

C. Hyper-skewness and hyper-kurtosis of net baryon-number fluctuations
The fifth-and sixth-order cumulants are related to the corresponding fifth-and sixth-order standardized moments, i.e., the hyper-skewness, S H , and hyper-kurtosis, κ H . We consider here the cumulant ratios for fifth-and sixth-order cumulants of net baryon-number fluctuations, Unlike the ratios for skewness and kurtosis cumulants, the corresponding ratios involving fifth-and sixth-order cumulants are negative already at μ B ¼ 0 in a broad temperature interval in the vicinity of T pc ð0Þ and become smaller with increasing μ B . This reflects the properties of the sixth-and eighth-order cumulants shown in Fig. 1.
The μ B dependence of the cumulant ratios R B 51 and R B 62 follows a pattern similar to that of the skewness and kurtosis ratios. In particular, at LO both ratios are almost identical and the NLO correction to R B 62 is about a factor of 3 larger than that for R B 51 . Like in the case of the corresponding relations for the skewness and kurtosis ratios these relations simply result from the structure of Taylor expansions for odd and even cumulants [19]. The relations are exact for expansions at vanishing μ Q and μ S and apparently they are not much altered in the strangeness neutral case n S ¼ 0 with n Q =n B ¼ 0.4. A fit to the N τ ¼ 8 lattice QCD results for the difference R B 62 − R B 51 at μ B ¼ 0 yields 0.029 (9). While statistical errors are strongly correlated between the fifth-and sixth-order cumulants they are large for each of these cumulants individually. For this reason we only present results for these cumulants obtained on lattices with temporal extent N τ ¼ 8 and evaluate the NLO corrections only for μ B =T ≤ 0.8. NLO results for R B 51 ðT; μ B Þ and R B 62 ðT; μ B Þ are shown in Fig. 8. Obviously NLO corrections for these ratios are negative and substantially larger than those in the skewness and kurtosis ratios. In the vicinity of the pseudocritical temperature the difference between LO and NLO results at μ B =T ¼ 0.8 is about an order of magnitude larger in R B 51 ðT; μ B Þ than in R B 31 ðT; μ B Þ. This is also the case when comparing R B 62 ðT; μ B Þ with R B 42 ðT; μ B Þ. The magnitude and sign of the NLO corrections to fifthand sixth-order cumulants in relation to corresponding results for the third-and fourth-order cumulants is evident from the structure of the corresponding Taylor expansion coefficients. It is easy to see this in Taylor expansions performed at μ Q ¼ μ S ¼ 0. In this case one has, for instance, As can be deduced from Fig. 1, despite the large errors on current results for χ B 8 , the cumulants χ B 6 and χ B 8 are both negative in the vicinity of the pseudocritical temperature; however the absolute value of the eighth-order cumulant is about an order of magnitude larger. This results in the much larger NLO correction to the expansion of χ B 6 ðT; μ B Þ. Although the expansions of all cumulants χ B n ðT; μ B Þ will have the same radius of convergence it is apparent that expansions for higher-order cumulants will converge more slowly. Higher-order corrections to χ B 5 ðT; μ B Þ and χ B 6 ðT; μ B Þ will thus be needed to arrive at firm conclusions on the behavior of these cumulants close to μ B =T ≃ 1. For μ B =T ≃ 0.3; however, the NLO correction is about an order of magnitude smaller and thus of similar magnitude as the NNLO correction to χ B 3 ðT; μ B Þ and χ B 4 ðT; μ B Þ at μ B =T ≃ 1. For small values of the baryon chemical potential and μ S ¼ μ Q ¼ 0 we thus may extend the result on the ordering of cumulant ratios stated in Eq. (3) and also include results for the fifth-and sixth-order cumulant ratios,

IV. BARYON-NUMBER FLUCTUATIONS ON THE PSEUDOCRITICAL LINE AND THE CUMULANTS OF NET PROTON-NUMBER FLUCTUATIONS
In this section we compare results on higher-order cumulants of net proton-number fluctuations, obtained by the STAR Collaboration during BES-I at RHIC [18,20], with our results for cumulants of net baryon-number fluctuations calculated in QCD on the pseudocritical line given in Eq. (1). The pseudocritical line shows only a rather weak dependence on μ B . The Oðμ 4 B Þ correction to T pc ðμ B Þ is found to be zero within errors [19]. For μ B ≤ T pc ð0Þ it changes from T ¼ 156.5ð1.5Þ to 154.5(2.0) MeV. This range of temperatures is well covered by the results for cumulant ratios as a function of μ B evaluated at fixed values of the temperature that have been shown in the previous section.
In Fig. 9 we show results for R B 12 ðT pc ðμ B Þ; μ B Þ on the pseudocritical line and compare with results obtained from noninteracting HRG model calculations that utilize hadron resonance gas spectra as listed in the Particle Data Tables [25] as well as spectra calculated in quark models [26,27]. As can be seen in Fig. 9 HRG model calculations for R B 12 agree well with QCD results obtained on the pseudocritical line up to about μ B =T ≃ 0.8 or μ B ≃ 125 MeV. This suggests that the use of low-order HRG cumulants, in particular the mean of hadron distributions (hadron yields) that are used experimentally to determine freeze-out parameters, may be appropriate at small values of the baryon chemical potential or small net baryon-number densities. The HRG model estimates of freeze-out parameters [12] suggest that the range of baryon chemical potentials μ B =T ≲ 1 corresponds to thermal conditions at freeze-out generated in heavy-ion experiments at beam energies ffiffiffiffiffiffiffi ffi s NN p ≳ 27 GeV. Figure 9 suggests that below this value of ffiffiffiffiffiffiffi ffi s NN p HRG model determinations of baryon chemical potentials differ from QCD determinations by more than 10%. It thus may be useful to eliminate μ B in favor of a directly accessible physical observable, e.g., R B 12 . At least for μ B ≲ 200 MeV truncation errors in the Taylor expansion of the first two cumulants, the mean and variance, as well as lattice discretization errors are small. The continuum limit extrapolation for R B 12 ðT pc ðμ B Þ; μ B Þ, shown in Fig. 9 thus does not suffer from truncation errors in the Taylor  series at least up to μ B =T ¼ 1.2. It is a monotonically rising function 3 of μ B . This allows to replace the chemical potential in an analysis of higher-order cumulant ratios in favor of R B 12 . We have done so for the comparison of higher-order cumulant ratios calculated in lattice QCD on the pseudocritical line with experimental data on cumulants of net proton-number fluctuations. In Fig. 10 we show the skewness and kurtosis ratios, R B 31 and R B 42 , on the pseudocritical line as a function of R B 12 , which also has been evaluated on the pseudocritical line. Similar results for the hyper-skewness and hyper-kurtosis ratios are shown in Fig. 11.
In Fig. 10 we show lattice QCD results up to R B 12 ¼ 0.75, which corresponds to μ B ¼ T pc ðμ B Þ ≃ 154.5 MeV. The width of the bands shown in the figure reflect the error on T pc ðμ B Þ as given in Eq. (1) as well as the error on the NNLO and continuum limit estimates for R B 31 and R B 42 . Note that the upper ends of these error bands correspond to the lower temperature, i.e., T ¼ 155 MeV at μ B ¼ 0 and T ≃ 152.5 MeV at μ B =T ¼ 1.
Also shown in this figure are results for the skewness and kurtosis ratios of net proton-number fluctuations obtained by the STAR Collaboration [18,20]. These ratios are plotted versus the measured ratio of the mean and variance of net proton-number fluctuations, which is taken as a proxy for the net baryon-number cumulant ratio 4 R B 12 .
As the experimentally determined skewness ratio of net proton-number fluctuations has a rather weak dependence on R P 12 and also the QCD result for R B 31 has a weak dependence on R B 12 , it obviously is not of much importance for the comparison of data and lattice QCD calculations whether R P 12 equals R B 12 or only is a proxy within say 10-20%. More relevant is the question to what extent the magnitude of R P 31 is a good approximation 5 for R B 31 . A direct comparison between R P 31 and R B 31 , as shown in Fig. 10, suggests that freeze-out happens in the vicinity of but below the pseudocritical temperature. In fact, as can be seen in Figs. 4 and 7, the ratios R B 31 and R B 42 are decreasing functions of the temperature. Experimental data for R P 31 lying above the theoretical band for R B 31 , evaluated on the pseudocritical line, thus suggest that freeze-out happens at a lower temperature.
Although errors on experimental results for the kurtosis ratio R P 42 are large, they are thermodynamically consistent with the data on the skewness ratio as pointed out already in our earlier analysis [19]. This gets further support through recent high-statistics 6 data obtained by the STAR Collaboration at ffiffiffiffiffiffiffi ffi s NN p ¼ 54.4 GeV [20]. These data are shown in Fig. 10 [20]. 3 Note that this will no longer be the case when one comes close to a critical point, where χ B 2 is expected to diverge and R B 12 ðT pc ðμ B Þ; μ B Þ thus would approach zero. 4 In a noninteracting HRG with vanishing strangeness and electric-charge chemical potential the ratios of the mean and variance of net proton-number fluctuations and net baryonnumber fluctuations are identical. In the case of a strangeness neutral (n S ¼ 0 with n Q =n B ¼ 0.4), noninteracting HRG, however, the latter is about 10% smaller. 5 Many caveats for a direct comparison between net baryonnumber fluctuations calculated in equilibrium thermodynamics and net proton-number fluctuations measured in heavy-ion collisions have been discussed in the literature [10,13]. The lattice QCD results shown in Fig. 10 thus may be considered only as a starting point for a more refined analysis of the experimental data that may take into account effects arising from experimental acceptance cuts, the small size of the hot and dense medium, nonequilibrium effects etc. 6 The statistics at ffiffiffiffiffiffiffi ffi Also shown in Fig. 10 with dashed lines is a joint fit to the experimental data on R P 31 and R P 42 [18] for ffiffiffiffiffiffiffi ffi s NN p ≥ 19.6 GeV using a quadratic ansatz, already used in Ref. [19], with K 0 ≡ S 0 . Including the new data at ffiffiffiffiffiffiffi ffi s NN p ¼ 54.4 GeV yields a fit, consistent with Ref. [19], but further constrains the parameters. One finds S 0 ≡ K 0 ¼ 0.761ð20Þ, S 2 ¼ −0.077ð70Þ, and K 2 ¼ −0.54ð22Þ. From the continuum estimates of R B 31 and R B 42 at μ B ¼ 0 shown in Fig. 6 one finds that the value of S 0 corresponds to a freeze-out temperature of 153.5(2.0) MeV. This temperature range is consistent with an earlier determination of the freeze-out temperature that was based on a comparison of the mean-tovariance ratio of net electric-charge and net proton-number ratios obtained by the STAR and PHENIX collaborations [29,30] with corresponding lattice QCD calculations for net electric-charge and net baryon-number cumulant ratios [31]. We also note that the ratio of the curvature of R B 42 and R B 31 on the pseudocritical (freeze-out) line tends to be larger than 3, which also has been noted in our previous analysis of the skewness and kurtosis ratios [19].
While the experimental data on the skewness and kurtosis cumulant ratios of net proton-number fluctuations, obtained at ffiffiffiffiffiffiffi ffi s NN p ≥ 27 GeV, are consistent with results on net baryon-number cumulants calculated within equilibrium QCD thermodynamics, this is not the case for the preliminary data on sixth-order cumulants presented by the STAR Collaboration [20]. The still preliminary data at ffiffiffiffiffiffiffi ffi s NN p ¼ 200 and 54.4 GeV, taken from the 0-40% centrality class, are shown in Fig. 11 together with the NLO lattice QCD calculations. At both values of ffiffiffiffiffiffiffi ffi s NN p deviations from the NLO lattice QCD results are large and of similar magnitude. While it is conceivable that the NLO results at R B 12 ≃ 0.5 (or μ B =T ≃ 0.6) will receive sizable corrections at NNLO, this is not the case at R B 12 ≃ 0.15 (or μ B =T ≃ 0.3). It thus seems impossible to describe both data points within QCD equilibrium thermodynamics. We also note that a large positive χ B 10 is needed, for such a contribution to render the hyperkurtosis ratio positive at ffiffiffiffiffiffiffi ffi s NN p ¼ 54. 4 GeV.
As pointed out in the previous section the NLO corrections for the hyper-skewness ratio R B 51 are a factor of 3 smaller than those for the hyper-kurtosis ratio R B 62 . Truncation errors for the former series are thus expected to be less severe. Furthermore, this ratio will also be easier to determine experimentally with smaller statistical errors. It thus would be an important check on the thermodynamic consistency of higher-order cumulants to compare experimental data on R P 51 at ffiffiffiffiffiffiffi ffi s NN p ≥ 54.4 GeV with the NLO lattice QCD calculations presented here.

V. SUMMARY AND CONCLUSIONS
We have presented new results on the μ B dependence of up to sixth-order cumulants using our latest results on up to eighth-order cumulants calculated at vanishing chemical potentials. Using simulation results obtained on lattices of size 32 3 × 8 and 48 3 × 12 we further presented continuum limit estimates for up to fourth-order cumulant ratios. For this analysis we used results from NNLO expansions of cumulants in the baryon chemical potential for strangeness neutral systems, n S ¼ 0 at an electric-charge to baryonnumber ratio n Q =n B ¼ 0.4. Systematic effects arising from the truncation of Taylor series for the skewness and kurtosis ratios were shown to be small for μ B =T ≤ 1, i.e., for the range of chemical potentials that can be probed in heavy-ion collisions in a range of beam energies ffiffiffiffiffiffiffi ffi s NN p ≥ 27 GeV. A comparison of the results on ratios of up to fourth-order cumulants of net baryon-number fluctuations calculated in equilibrium QCD thermodynamics with corresponding cumulants of net proton-number fluctuations yields quite good agreement. This suggests that the latter are consistent with reflecting the imprint of thermodynamical fluctuations at a temperature close to but below the pseudocritical temperatures T pc ðμ B Þ. The particularly good agreement between lattice QCD calculations and the high-statistics experimental data for up to fourth-order cumulants at ffiffiffiffiffiffiffi ffi s NN p ≥ 54.4 GeV suggests that this conclusion could be further strengthened, if data with similarly high statistics also becomes available at other beam energies in the range ffiffiffiffiffiffiffi ffi s NN p ≥ 27 GeV. We also presented first results from a NLO calculation of fifth-and sixth-order cumulants and showed that the hyperskewness and hyper-kurtosis ratios R B 51 and R B 62 are negative at low values of μ B =T and temperatures in the vicinity of T pc ðμ B Þ. This is at odds with preliminary data obtained by the STAR Collaboration at ffiffiffiffiffiffiffi ffi s NN p ≥ 54.4 GeV for the sixth-order cumulant ratio, R P 62 , of net protonnumber fluctuations, which was found to be positive and close to unity. However, on the one hand corrections to the LO result for these cumulants, calculated in lattice QCD, are large already at μ B ≃ 0.5. This makes a calculation of NNLO corrections for these cumulants desirable. On the other hand, the experimental determination of sixth-order cumulant ratios is known to require high statistics and current experimental data may be statistics limited. We pointed out that a measurement of ratios of fifth-and firstorder cumulants would be very helpful as this ratio can be better controlled experimentally and suffers less from truncation effects in NLO lattice QCD calculations.
All data from our calculations, presented in the figures of this paper, can be found at https://pub.uni-bielefeld.de/ record/2941824 [32]. (iii) the INCITE program at Argonne Leadership Computing Facility, a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-06CH11357; and (iv) the USQCD resources at the Thomas Jefferson National Accelerator Facility. This research also used computing resources made available through: (i) a PRACE grant at CINECA, Italy; (ii) the Gauss Center at NIC-Jülich, Germany; and (iii) the Nuclear Science Computing Center at Central China Normal University and (iv) the GPU-cluster at Bielefeld University, Germany.