Higher order fluctuations and correlations of conserved charges from lattice QCD

We calculate several diagonal and non-diagonal fluctuations of conserved charges in a system of 2+1+1 quark flavors with physical masses, on a lattice with size $48^3\times12$. Higher order fluctuations at $\mu_B=0$ are obtained as derivatives of the lower order ones, simulated at imaginary chemical potential. From these correlations and fluctuations we construct ratios of net-baryon number cumulants as functions of temperature and chemical potential, which satisfy the experimental conditions of strangeness neutrality and proton/baryon ratio. Our results qualitatively explain the behavior of the measured cumulant ratios by the STAR collaboration.


Introduction
One of the most challenging goals in the study of Quantum Chromodynamics (QCD) is a precise mapping of the phase diagram of strongly interacting matter. First principle, lattice QCD simulations predict that the transition from hadrons to deconfined quarks and gluons is a smooth crossover [1][2][3][4][5][6], taking place in the temperature range T 145 − 165 MeV. Lattice simulations cannot presently be performed at finite density due to the sign problem, thus leading to the fact that the QCD phase diagram is still vastly unexplored when the asymmetry between matter and antimatter becomes large.
With the advent of the second Beam Energy Scan (BES-II) at the Relativistic Heavy Ion Collider (RHIC), scheduled for 2019-2020, there is a renewed interest in the heavy ion community towards the phases of QCD at moderate-to-large densities. A rich theoretical effort is being developed in support of the experimental program; several observables are being calculated, in order to constrain the existence and location of the QCD critical point and to observe it experimentally.
Fluctuations of conserved charges (electric charge Q, baryon number B and strangeness S) are among the most relevant observables for the finite-density program for several reasons. One possible way to extend lattice results to finite density is to perform Taylor expansions of the thermodynamic observables around chemical potential µ B = 0 [7][8][9][10][11]. Fluctuations of conserved charges are directly related to the Taylor expansion coefficients of such observables, thus, they are needed to extend first principle approaches to the regions of the phase diagram relevant to RHIC. An other popular method to extend observables to finite density is the analytical continuation from imaginary chemical potentials [12][13][14][15][16]. The agreement between the analytical continuation and Taylor expansion was shown for the transition temperature by Bonati et al in Ref. [17].
Fluctuations can also be measured directly, and a comparison between theoretical and experimental results allows to extract the chemical freeze-out temperature T f and chemical potential µ Bf as functions of the collision energy [18][19][20][21][22]. Such fluctuations have been recently calculated and extrapolated using the Taylor method in Ref. [23]. Finally, higher order fluctuations of conserved charges are proportional to powers of the correlation length and are expected to diverge at the critical point, thus providing an important signature for its experimental detection [9,24,25].
In this paper, we calculate several diagonal and non-diagonal fluctuations of conserved charges up to sixth-order and give estimates for higher orders, in the temperature range 135 MeV ≤ T ≤ 220 MeV, for a system of 2+1+1 dynamical quarks with physical masses and lattice size 48 3 × 12. We simulate the lower-order fluctuations at imaginary chemical potential and extract the higher order fluctuations as derivatives of the lower order ones at µ B = 0. This method has been successfully used in the past and proved to lead to a more precise determination of the higher order fluctuations, compared to their direct calculation [26,27]. The direct method (see e.g. [7]) requires the evaluation of several terms and is affected by a signal-to-noise ratio which is decreasing as a power law of the spatial volume V , with an exponent that grows with the order of the susceptibility.
We also construct combinations of these diagonal and non-diagonal fluctuations in order to study the ratio of the cumulants of the net-baryon number distribution as functions of temperature and chemical potential by means of their Taylor expansion in powers of µ B /T . We discuss their qualitative comparison with the experimental results from the STAR collaboration, as well as the validity of the truncation of the Taylor series.
The paper is organized as follows: we first discuss the use of imaginary chemical potentials in Section 2. Section 3 gives details on the lattice setup, on the fitting procedure, on its generalization for cross-correlators, and finally on the error estimation. The phenomenological results for the ratios of kurtosis, skewness and variance of the baryon number are presented in Section 4. Conclusions and outlook are discussed in Section 5, while in the Appendix we present all diagonal and non-diagonal fluctuations needed to construct the cumulant ratios shown in Section 4, and give additional technical details. are given by The observables we are looking at are the derivatives of the free energy with respect to the chemical potentials. Since the free energy is proportional to the pressure, we can write: These are the generalized fluctuations we calculated around µ = 0 in our previous work [28]. The fermion determinant det M (µ) is complex for real chemical potentials, prohibiting the use of traditional simulation algorithms. For imaginary µ, however, the determinant stays real. The chemical potential is introduced through weighted temporal links in the staggered formalism: Thus, an imaginary µ translates into a phase factor for the antiperiodic boundary condition in the Dirac operator. Due to the Z(3) symmetry of the gauge sector, there is a non-trivial periodicity in the imaginary quark chemical potential µ q → µ q + i(2π/3)T , which translates to the baryochemical potential as µ B → µ B + i2πT , the Roberge-Weiss symmetry. This is independent of the charge conjugation symmetry µ B ↔ −µ B . As a result, e.g. for the imaginary part of the baryon density: At µ B = iπT there is a first order phase transition at all temperatures above the Roberge-Weiss critical end point T RW [29]. When µ B crosses iπT in the imaginary direction, the imaginary baryon density is discontinuous. This behaviour is illustrated in Fig. 1, where the imaginary baryon density as a function of the imaginary chemical potential is shown. At low temperature the Hadron Resonance Gas model predicts B ∼ sinh(µ B /T ), thus for imaginary values we expect a sine function below T c : Im B ∼ sin(Imµ B /T ). At temperatures slightly above T c , we observe that further Fourier components appear in addition to sin(Imµ B /T ) with alternating coefficients, these are consistent with a repulsive interaction between baryons [30]. At very high temperatures, on the other hand, B is a polynomial of µ B since the diagrams contributing to its ∼ µ 5 B and higher order components are suppressed by asymptotic freedom [31,32]. The Stefan-Boltzmann limit is non-vanishing only for two Taylor coefficients of Im B , giving Im B | µ B /T =iπ− = 8π/27. At finite temperatures above T RW this expectation value is smaller but positive. By Eq. (2.5), it implies a first order transition at µ B = iπT .
The order of the transition at T RW heavily depends on the quark masses [33,34]. For physical quark masses one obtains T RW = 208(5) MeV, and the scaling around the endpoint is consistent with the Ising exponents [35]. This implies that, for physical parameters, the transition is limited to µ B = iπT without any other structures between the imaginary interval [0, iπ) [33].
Thus, we have only the range µ/T ∈ [0, iπ) to explore the µ-dependence of the observables. Recent simulations in this range include the determination of the transition line, where the slope was determined on the negative side of the T − µ 2 B phase diagram. Using analyticity arguments, this coefficient gives the curvature of the transition line on the real T − µ B phase diagram [36][37][38]. Apart from the transition temperature, we used imaginary chemical potentials also to extrapolate the equation of state to real µ B [26], which serves as an alternative approach to the Taylor extrapolation [39]. In an recent study D'Elia et al. have used the low order fluctuations at imaginary chemical potentials to calculate generalized quark number susceptibilities [27].

Lattice setup
In this work we calculate high order fluctuations by studying the imaginary chemical potential dependence of various generalized quark number susceptibilities. We use a tree-level Symanzik improved gauge action, with four times stout smeared (ρ = 0.125) staggered fermions. We simulate 2 + 1 + 1 dynamical quarks, where the light flavors are tuned in a way to reproduce the physical pion and kaon masses and we set mc ms = 11.85 [40]. For the zero-temperature runs that we used for the determination of the bare masses and the coupling, the volumes satisfy Lm π > 4. The scale is determined via f π . More details on the scale setting and lattice setup can be found in [28].
Our lattice ensembles are generated at eighteen temperatures in the temperature range 135. . . 220 MeV. We simulate at eight different values of imaginary µ B given as: µ for j ∈ {0, 1, 2, 3, 4, 5, 6, 7}. In this work the analysis is done purely on a 48 3 × 12 lattice, we leave the continuum extrapolation for future work.
In terms of quark chemical potentials we generate ensembles with µ u = µ d = µ s = µ B /3. In each simulation point we calculate all derivatives in Eq. (2.2) up to fourth order. Thanks to our scan in Imμ B , we can calculate additional µ B derivatives. Ref. [27] uses various "trajectories" in the µ B − µ Q − µ S space, allowing the numerical determination of higher e.g. µ Q and µ S derivatives. We find relatively good signal for the µ Q and µ S derivatives by directly evaluating Eq. (2.2) within one simulation. We recently summarized the details of the direct calculation in Ref. [28].

Correlated fit with priors
We start with the analysis for χ B 2 (T ), χ B 4 (T ) and χ B 6 (T ). Our goal is to calculate these quantities at zero chemical potential, using the imaginary chemical potential data up to In this work we extract these derivatives at a fixed temperature. Results for different temperatures are obtained completely independently, an interpolation in temperature is not necessary at any point. Thus, the error bars in our results plot will be independent. The errors between the quantities χ B 2 (T ), χ B 4 (T ) and χ B 6 (T ) will be highly correlated, though, since these are extracted through the same set of ensembles at the given temperature. This correlation will be taken into account when combined quantities are calculated, or when an extrapolation to real chemical potential is undertaken.
Thus we consider the ensembles at a fixed temperature T . For each value of imaginary µ B = 0 we determine χ B 1 , χ B 2 , χ B 3 and χ B 4 from simulation, while for µ B = 0 only χ B 2 and χ B 4 can be used, since χ B 1 and χ B 3 are odd functions of µ B and therefore equal to zero. We make the ansatz for the pressure: where the Taylor expansion coefficients c n are related to the baryon number fluctuations χ B n by: n!c n = χ B n . Our data do not allow for an independent determination of c 8 and c 10 . Nevertheless, in order to have some control over these terms we make the assumption or in terms of the c n coefficients We can then rewrite our ansatz as where 1 and 2 are drawn randomly from a normal distribution with mean -1.25 and variance 2.75. We use the same distribution for all temperatures. In effect, our c 8 and c 10 coefficients are stochastic variables. There is sufficient probability to obtain coefficients that slightly break Eq. (3.3), should the data prefer larger than expected fluctuations, e.g. due to a nearby critical end point. The used distribution for 1,2 actually implements a prior for χ B 8 and χ B 10 . In the HRG model we know that , which is well represented by the prior distribution. At high temperatures χ 6 B and higher coefficients quickly approach zero, as obtained in Hard Thermal Loop results [41]. In the transition regime, the higher moments of the baryon fluctuations are dominated by the fact that the transition line is µ B -dependent. Starting from µ B = 0 at a fixed temperature between T c and T RW , a crossover line is developed as the imaginary chemical potential is introduced. The magnitude of higher order fluctuations in the transition regime can be estimated by a very simple observation. The behaviour of the quark density χ B , where κ is the curvature of the transition line in the µ B − T phase diagram. In this approximation, the only source of µ B -dependence is coming from the curvature of the transition line. Calculating the µ B derivatives gives a basic estimate for χ B 8 , which we used to tune the prior distribution.
For this ansatz we calculate the following derivatives, which are the actually simulated lattice observables: We perform a correlated fit for the four measured observables, thus obtaining the values of c 2 , c 4 and c 6 for each temperature, and the corresponding χ B 2 , χ B 4 and χ B 6 . We repeat the fit for 1000 random draws for 1 and 2 . The result is weighted using the Akaike Information Criterion [42]. Through these weights we get a posterior distribution from the prior distribution. Our final estimate for χ B 8 represents this posterior distribution. We do not show the posterior for χ B 10 , which is mostly noise. These results are shown in Fig. 2, together with an estimate of χ B 8 , related to χ B 4 by Eq. (3.4).

Cross-correlators
So far we only considered derivatives with respect to the baryonic chemical potential. In our previous, direct analysis in Ref. [28], the µ B -derivatives had larger errors than µ Q − or µ S −derivatives. For µ Q , the most noisy disconnected contributions come with smaller prefactors, while for µ S the disconnected contributions are small due to the heavier strange mass. Our approach was designed to improve the µ B -derivatives only. Therefore, the µ S and µ Q derivatives have to be simulated directly and without the support from the fit that we used in the µ B direction. Our result on χ QS jk improved only due to the increase in the statistics since [28].
On the other hand, baryon-strange and baryon-charge mixed derivatives do benefit from the imaginary µ B data. We simulate various χ B,Q,S i,j,k with the appropriate values of j and k and all possible values of i so that i + j + k ≤ 4. For each group of fluctuations with the same j and k we perform a fit analogous to the procedure described in Section 3.2.
Let's take the example of j = 1, k = 0. Our ansatz for cross-correlators is analogous to Eqs. (3.5)-(3.8): We truncated the expression at tenth order. The priors assume |χ BS 71 | |χ BS 31 | and |χ BS 91 | |χ BS 31 |, as it is certainly true at high temperature and within the HRG model. The prior distribution is wider than 1, we used the same mean and variance as in the channel with no µ S derivative. When we use Eq. (3.9) we take χ S 1 , χ BS 11 , χ BS 21 and χ BS 31 as correlated quartets for each imaginary chemical potential and determine the three free coefficients of Eq. (3.9). This fitting procedure is repeated 1000 times with random χ BS 71 /χ BS 31 and χ BS 91 /χ BS 31 coefficients. Again, using the Akaike weights we constrain the prior distribution. The resulting estimate for χ BS 71 along with the fit coefficients are shown in Fig. 3. The posterior for χ BS 91 is not only noisy, but it is probably heavily contaminated by the higher orders that we did not account for. The other channels with higher µ S or µ Q derivatives are obtained analogously. These are plotted in Appendix A.

Error Analysis
For a reliable comparison between experimental measurements and theoretical calculations, the error estimate is an important ingredient. Our statistical error is estimated through the jackknife method. For our systematic error there are several sources. We determine our systematic error by the histogram method described in [43], where each analysis is weighted with the Akaike information criteria. We include the influence of the number of points in the µ B direction, by either including or ignoring the data from our highest value of µ B . A very important source for our systematic error is the influence of the higher order contributions in µ B . This effect was estimated by adding the higher order terms with prefactors 1 and 2 as described in Section 3.2. We consider 1000 different pairs and add the different analyses to our histogram. The width of the histogram using Akaike weights corresponding to the fit quality gives the systematic errors for the fit coefficients, and from the same histogram we obtain the posterior distributions for 1 . The physical quantities that are constrained only by the posterior distribution are plotted with green symbols.
These histograms are built independently for each number (j and k) of µ S and µ Q derivatives. When calculating the systematics for the cumulant ratios (Section 4) we need to calculate different combinations of diagonal and non-diagonal fluctuations from the available analyses. Though these fits (corresponding to the same temperature) are carried out separately we keep track of the statistical correlation, by maintaining the jackknife ensembles throughout the analysis. The correct propagation of systematic errors is a more elaborate procedure. When χ BSQ ijk coefficients are combined with different j, k pairs, different histograms have to be combined. If we had only two variables to combine, each of the 2000 first fit variants should be combined with each of the 2000 second fit variants and use the product of the respective probability weights. Instead, we combine the fit results by drawing 'good' fits by importance sampling from each histogram independently. In this way, O(100) random combinations of χ BSQ ijk results already give convergence for each discussed quantity and its error bar. For the results in this paper we used 1000 such random combinations. This procedure assumes that between different j, k pairs the prior distribution is uncorrelated.

Phenomenology at finite chemical potential
For a comparison with heavy ion collision experiments, the cumulants of the net-baryon distribution are very useful observables. The first four cumulants are the mean M B , the variance σ 2 B , the skewness S B and the kurtosis κ B . By forming appropriate ratios, we can cancel out explicit volume factors. However, the measured distributions themselves may still depend on the volume, which one should take into account when comparing to experiments.
Heavy ion collisions involving lead or gold atoms at µ B > 0 correspond to the following situation n S = 0 n Q = 0.4 n B . (4.1) For each T and µ B pair, we have to first calculate µ Q and µ S that satisfy this condition. The resulting µ Q (µ B ) and µ S (µ B ) functions, too, can be Taylor expandend [19,20], introducing We investigate three different ratios of cumulants: The µ B -dependence of the χ B i (T,μ B ) can again be written as a Taylor series: The χ coefficients that we determined in Section 3 include derivatives up to sixth order, and we have estimates for the eighth order, too. The fit coefficients corresponding to the tenth order are likely to be contaminated by higher orders, that we did not include into the ansatz. These χ BQS ijk coefficients, however, are given for j + k ≤ 4, which is the highest order that we used in µ Q and µ S . This list of coefficients allows us to calculate the r B,k ij coefficients from Equations (4.4), (4.5) and (4.6). The results for the r B,k ij coefficients are shown in Figures 4, 5 and 6. We confirm the observation from Ref. [23] that the coefficient r B,2 42 has a similar temperature dependence as r B,2 31 but it is ∼ 3 times larger in magnitude. For higher order coefficients, higher order derivatives in µ S and µ Q are needed. The direct simulations have a rapidly increasing error with the order of the derivative, and very large statistics would be needed to improve our calculations at this point. Another possibility would be to simulate new ensembles with finite µ S and µ Q and do a similar fit as for the µ B direction. This approach has been used in [27].
After calculating the Taylor coefficients for S B σ 3 B /M B and κ B σ 2 B , we use these results to extrapolate these quantities to finite chemical potential. They are shown in Figure 7. In the left panel, S B σ 3 B /M B is shown as a function of the chemical potential for different temperatures. The Taylor expansion for this quantity is truncated at O(μ 2 B ). The black points in the figure are the experimental results from the STAR collaboration from an analysis of cumulant ratios measured at mid-rapidity, |y| ≤ 0.5, including protons and antiprotons with transverse momenta 0.4 GeV ≤ p t ≤ 2.0 GeV [44,45]. The beam energies were translated to chemical potentials using the fitted formula of Ref. [46]. Even if we do not quantitatively compare the lattice bands to the measurements to extract the freeze-out parameters, as experimental higher order fluctuations might be affected by several effects of non-thermal origin and our lattice results are not continuum extrapolated, we notice that the trend of the data with increasing µ B can be understood in terms of our Taylor expansion.
In the STAR collaboration with transverse momentum cut 0.4 GeV≤ p t ≤ 2.0 GeV [44,45]. Notice that, due to the fact that the r B,4 42 is positive in the range 160 MeV≤ T ≤195 MeV, we observe a non-monotonic behavior in κ B σ 2 B for T = 160 MeV at large chemical potentials. By comparing the two different truncations of the Taylor series we can conclude that, as we increase the temperature, the range of applicability of our Taylor series decreases: while at T = 150 MeV the two orders agree in the whole µ B /T range shown in the figure, at T = 160 MeV the central line of the next-to-next-to-leading order bends upwards and is not contained in the next-to-leading order band. To make the NLO prediction precise substantially more computer time would be needed.

Conclusions and outlook
In this manuscript, we have calculated several diagonal and non-diagonal fluctuations of electric charge, baryon number and strangeness up to sixth-order, in a system of 2+1+1 quark flavors with physical quark masses, on a lattice with size 48 3 × 12. The analysis has been performed simulating the lower order fluctuations at zero and imaginary chemical potential µ B , and extracting the higher order fluctuations as derivatives of the lower order ones at µ B = 0. The chemical potentials for electric charge and strangeness have both been set to zero in the simulations. From these fluctuations, we have constructed ratios of baryon number cumulants as functions of T and µ B , by means of a Taylor series which takes into account the experimental constraints n S = 0 and n Q = 0.4 n B . These ratios qualitatively explain the behavior observed in the experimental measurements by the STAR collaboration as functions of the collision energy.
We focused on observables (baryon distribution, ratios of cumulants) that are less sensitive to lattice artefacts. An obvious extension of our work will be the use of finer lattices and a continuum extrapolation. The other extension is to use a two-or even three-dimensional mapping of the space of the imaginary chemical potentials using nonvanishing µ S and µ Q . That would not only improve the µ S − and µ Q −derivatives, but would allow us to study the melting of states with various strangeness and electric charge quantum numbers. Our first study in this direction using strangeness chemical potentials was published in Ref. [47].

A Results for the correlators
In this Appendix we present the non-diagonal fluctuations of conserved charges needed to construct the cumulant ratios at finite chemical potential µ B , satisfying the constraints n S = 0 and n Q = 0.4 n B .
Like we did for the diagonal χ B i , we simulate lower order fluctuations at finite imaginary chemical potential and extract the higher order fluctuations as derivatives of the lower order ones at µ B = 0: in particular, we simulate various χ B,Q,S i,j,k with the appropriate values of j and k and all possible values for i so that and extract the corresponding χ B,Q,S i,j,k with i + j + k ≤ 6 and an estimate for i + j + k = 8 and sometimes even i + j + k = 10. By estimate (shown in green) we mean the posterior distribution that we get for the two highest orders when using priors, as discussed in the main text. In total we need 15 channels to obtain all the necessary terms.
In the following plots we show these results organized by the number of charge derivatives (j) in Figs. [8][9][10][11][12]. It is notoriously difficult to calculate charge correlators using staggered fermions [28]. Correlators that are not protected by a baryon derivative are affected by significant discretization errors. It is understood in the HRG model context that discretization errors mostly affect the contributions from pions and kaons. Staggered lattice effects introduce the highest relative errors for the lightest mesons. Luckily, however, quantities with such discretization effects come with a small pre-factor into the final formulas of Eqs.˙(4.4)-(4.6). If we had a complete isospin symmetry (factor 0.5 between n Q and n B in Eq. (4.1)) then electric charge correlators would play no role at all in the extrapolation of baryon fluctuations.

B Statistics and lattice details
In Table 1 we give the number of analyzed configurations per ensemble. The simulation parameters and the details of the analysis are given in Ref. [28].
The determination of the µ derivatives follows the lines of Ref. [7,28]. We calculate four quantities per configuration and per quark mass Here M j is the fermion matrix corresponding to the j-th quark mass in the system. M and M indicate the first and higher order derivatives with respect to the quark chemical potential. For this simple staggered action higher order derivatives are equal to lower order ones, M = M and M = M by construction. These traces are calculated using the standard stochastic method, by calculating the effect of the matrices on random sources. At finite (imaginary) chemical potentials we used 4 × 256 Gaussian random sources for the light quarks and 4 × 128 sources for the strange quarks. The analysis was accelerated by calculating 256 eigenvectors of the Dirac operator first. These eigenvectors were then fed into an Eig-CG algorithm. Using the isospin symmetry (m u = m d ), the ABCD traces can be used to calculate the χ uds derivatives with the following formulas: . Results containing one electric charge derivative on the various correlators on our 48 3 ×12 lattice as functions of the temperature. Green data points denote our estimates for the high orders, these were fitted using a prior distribution. The red curves are the HRG model results.  Figure 10. Results containing two electric charge derivatives on the various correlators on our 48 3 × 12 lattice as functions of the temperature. Green data points denote our estimates for the high orders, these were fitted using a prior distribution. The red curves are the HRG model results. Charge correlators without baryon derivative (here χ QS 22 ) are expected to have significant discretization errors.
(B.14)  Figure 12. χ Q 4 , χ BQ 24 and estimate for χ BQ 44 as functions of the temperature. The quantity χ Q 4 has severe cut-off effects on this lattice [28]. The red curves are the HRG model results.
If the listed products of the A, B, C, D traces are calculated as products of the stochastic estimators, a bias could be introduced. Thus, in products different random vectors have to be used in each factor. Alternatively, the expectation value of the bias has to be subtracted. The last step is to express the derivatives in terms of µ B , µ Q and µ S in Eq. (2.2) using Eqs. (2.1), which is a straightforward exercise.  Table 1. Statistics of our simulations on the 48 3 × 12 lattice. We list the number of stored and analyzed gauge configurations. These configurations were separated by ten Rational Hybrid Monte Carlo updates.