High order quark number susceptibilities in hot QCD from lattice EQCD

Building on the experience of [1], we develop a formalism to construct operators for higher derivatives of the pressure in hot QCD with respect to the quark chemical potential $\mu$. We provide formulae for the operators up to the sixth derivative, and obtain continuum-extrapolated results from lattice EQCD at zero and finite $\mu$ and at six different pairs of temperature T and number of massless quark flavors $n_f$. Our data is benchmarked against full-QCD lattice and perturbative results, allowing to judge the quality of the perturbative series expansion in EQCD and the dimensional reduction procedure as a whole.


Introduction
It has been discovered at RHIC [2,3] and LHC [4][5][6][7] that the Quark-Gluon Plasma (QGP) behaves effectively as an almost perfect fluid. Therefore, it seems natural to attempt modeling the QGP's dynamics with hydrodynamic equations [8,9]. Beyond equations of continuity for all conserved charges and an equation of motion for the velocity field, a thermodynamic equation of state (EOS) is required. While experimental measurements of the equation of state are underway, our work will focus on the theoretical side. A broad variety of methods has been applied to study the equation of state of the QGP, i.e. the pressure p as a function of the temperature T and the quark chemical potential µ, for instance holographic methods [10], functional methods [11,12], perturbation theory [13,14], and lattice simulations [15,16]. Our approach will concentrate on the two latter, methodologically more conservative approaches.
Perturbation theory relies on the system being weakly coupled, meaning that the temperature must be much higher than the renormalization scale Λ. In particular, the temperature must be far above the pseudocritical temperature T c of the predicted crossover between the hadronic phase of nuclear matter and the QGP phase at µ = 0. On the other hand, lattice computations tackle the path integral numerically. For convergence it is important that the measure of the path integral is strictly positive. A non-vanishing quark chemical potential, however, manifests in an imaginary part of the action [17]. The measure of the path integral becomes not strictly positive, and one runs into the famous sign problem of QCD which makes direct simulations at µ = 0 in practice impossible. Numerous ways have been developed to circumvent that problem, one of which is the Taylor expansion of the pressure in terms of the quark chemical potential µ around µ = 0 [18][19][20]. The Taylor coefficients itself are computed as increasingly complicated correlation functions at µ = 0. Currently, there are numerical results up to the sixth order in µ [21]. From the ratio of the coefficients, one can also estimate the radius of convergence of the Taylor series and therefore infer a lower bound on the location of the conjectured tricritical point on the µ-axis of the QCD phase diagram.
At temperatures of a few times the crossover temperature T c , QCD thermodynamics does not yet behave fully perturbatively. The reason for that are the gluon Matsubara zero modes which contribute a factor of the inverse strong coupling constant 1/g to each closed loop and therefore modify the g 2 -suppression per loop to g 2 /g = g [22]. Fortunately, these problematic modes can be separately treated in a three-dimensional effective theory, 'electrostatic QCD' (EQCD) [23][24][25]. Beyond the advantage of allowing for a separate treatment and resummation of the Matsubara zero modes, EQCD also has the upside of a much milder sign problem. In fact, an analytic continuation of the quark chemical potential is possible, enabling us to solve EQCD directly on the lattice at moderate non-vanishing µ.
The goal of this work is twofold: First, we compute Taylor coefficients from EQCD at vanishing chemical potential in order to compare to both lattice results at temperatures slightly above T c and perturbative predictions at very high temperature. The higher cumulants have proven to be an especially well-suited quantity to judge the lower end of the range of validity of EQCD in the past [1]. As even higher cumulants are expected to be even less dominated by ultraviolet effects, this should also hold for the higher cumulants. This part also serves as a crosscheck of our work. Second, we give numerical values of the quark number susceptibility up to sixth order directly simulated at nonvanishing chemical potential µ, a region that is still inaccessible to traditional lattice simulations.
The paper is organized as follows: We give the formulae for the dimensional reduction procedure from full QCD to EQCD in Sec. 2 and specify our scenarios of interest. A recipe for the derivation of operators for the higher cumulants and how to analytically continue them to make them accessible to lattice simulations is given in Sec. 3. Sec. 4 contains the details of our lattice implementation, in particular how the operators presented in the previous section translate to the lattice. Numerical results, for second to sixth derivatives with respect to the quark chemical potential µ are given in Sec. 5, and compared to predictions from perturbation theory and four-dimensional lattice simulations. We conclude with Sec. 6. Appendices contain a table with our simulation parameters, analytic continuation of the derivative expressions beyond the third order, and comprehensive tables with values for all derivatives at all temperatures.

Thermodynamics from Electrostatic Quantum Chromodynamics
At high temperatures thermal four-dimensional QCD reduces to electrostatic QCD, which only treats the gluonic Matsubara zero mode dynamically and hides all other modesgluonic and fermionic -in effective field theory parameters. The action of EQCD reads with the now dimensionful gauge coupling g 2 3d . Since all derivatives in time direction vanish, gauge freedom does not prevent the former A 0 field from acquiring a screening mass m 2 D . Therefore, the zero component of the four-dimensional gauge field A 0 turns into SU(3) adjoint representation scalar field Φ = φ a T a , where T a = λ a /2, a = 1 . . . 8 are the Gell-Mann matrices. The gluon self-interaction generates a four-point interaction proportional to λ. The term cubic in Φ explicitly violates the otherwise valid mirror symmetry Φ → −Φ. It can be directly linked to a nonvanishing quark chemical potential µ in four dimensions.
Due to its three dimensional nature, EQCD is a super-renormalizable theory, meaning that all correlation functions can be rendered finite at all orders with a finite number of counterterms. In particular, only the screening mass m 2 D in 2.1 receives counterterms (at 1 and 2 loop order), which therefore is the only scale-dependent parameter. Consequently, one can use the dimensionful gauge coupling g 2 3d to set the scale and proceed with the dimensionless ratios (2.2)

EFT parameter matching
The values of the EQCD parameters can be rigorously derived from a perturbative matching computation [26,27]. The starting point is the 1 loop running coupling of full QCD Consulting the overview over the state of the matching in Sec. 6.3.1 of [28], we find that the precision bottleneck is the z-matching, which is only available at O(g 3 ). Therefore, we can express the explicit dependence of the EFT parameters on g 2 through x, and truncate the expansion of the other parameters in powers of x at O(x), since parametrically x ∼ g 4 /g 2 = g 2 .
Hereμ i = µ i /T , and µ i is the quark chemical potential for flavor i. Relevant scales for the matching areμ and For the (perturbative) scale-setting, we use Λ MS = 341 MeV determined in [29] for threeflavor QCD. We investigate EQCD corresponding to the six scenarios of full QCD in Tab. 1. An analytic continuation iz → z will finally allow us to get rid of the prefactor of i in front of the cubic term in (2.1). This procedure was found to be highly superior to the more naive continuation ofμ [1], since we only continue the linear contributions inμ, which appear in z, and leave the quadratic, i.e. real, contributions byμ to y (see (2.5)) untouched. For the time being, especially for the derivation of the investigated operators, we keep the cubic term in its original complex form, and perform the analytic continuation of action and operators just before discretizing the theory on the lattice. For each of the scenarios in Table 1

Matching the pressure to full QCD
Our goal is to study different susceptibilities, i.e. derivatives of the pressure, in full QCD using simulations of EQCD. To this end, we do not only have to match the parameters in the EQCD action to full QCD, but also the pressure itself. The philosophy for this is quite simple; we take the perturbative pressure of QCD, subtract its perturbative EQCD counterpart and substitute it with the lattice EQCD result. The remaining perturbative part is UV-dominated and can be reliably evaluated in perturbation theory: Depending on which derivative we want to investigate, we take derivatives of the above expression with respect to a given number of µ i . The four-dimensional pressure was computed perturbatively up to O(g 5 ) in ref. [14]. This computation also contains the threedimensional pressure, which is divergent by itself, though. These divergences cancel against divergences in the hard sector in the four-dimensional case. In the purely three-dimensional case, however, we need to introduce a divergent cosmological constant into the EQCD action that cures these divergences.

Higher-order cumulants in EQCD
Beyond the second order derivatives studied in [1], one can also investigate higher derivatives of the pressure with respect to different quark chemical potentials, sometimes also called cumulants. Since we only consider massless quarks, they are degenerate and we do not need to specify the name of the quark, e.g.
and conversely for higher derivatives. Even though a brute force derivation of the corresponding correlation functions in EQCD is possible, it can be dramatically simplified with a little additional work.

Derivation of higher order cumulants
We calculate derivatives of the pressure, i.e.
with the free energy density f and the volume V . Recalling that only fully connected correlation functions contribute to the free energy -and consequently the pressure -helps simplifying the calculation drastically. Subtractions of partly connected correlation functions can therefore be dropped in the following, and re-instated in the end. This also means that we only have to keep track of all derivatives acting directly on Z, since the derivatives acting on the 1/Z-normalization of correlation functions generate the partly connected subtractions. Moreover, we need to take a closer look at the (μ-dependent) part of the EQCD action Derivatives can only act on Z in two different ways: either as a single derivative on the exponential function or as a double derivative 1 Note the minus sign due to the weight of the exponential in the Euclidean path-integral being exp(−S EQCD ). Calculating ever-higher derivatives of the pressure thus amounts to combining the two above expressions in all possible ways maintaining the symmetry in the indices. Similar to the subtractions, it is sufficient to keep only specific combinations of the indices and re-instate the symmetry of the derivative in the end, bearing in mind that derivatives acting on continuous functions commute. Since the dimensional reduction to EQCD only works away from any phase transition, i.e. a non-analyticity in the free energy y, anyway, it is safe to assume that derivatives with respect to the quark chemical potentialμ acting on the pressure commute. As a first pedagogical example, we re-derive the quadratic quark number susceptibility, as used in [1]. Since it is a second derivative, (3.4) can only contribute as and (3.5) as where the symmetrization in the indices i and j is already intact and does not yet need to be enforced at this stage. We combine the two above terms to with the help of the condensates Let us caution the reader once more that the subtractions of the partly connected contributions have to be re-instated in the end. Keeping track of all the subtractions throughout the derivation does not provide any further insights about the underlying mechanisms. As a first non-trivial example, we present the derivation of the expression for the cubic derivative of the pressure with respect to the chemical potential. Once more, this can only receive two possible contributions: where the dots will denote the symmetrizations in the respective indices in the following. We can combine these two possibilities to where we defined the condensates Considering (3.11), we observe that by settingμ = 0, only K 2,2 and K 3,4 survive. These two condensates, in turn, are odd in Φ 3 , therefore they should be consistent with 0 numerically at vanishing z. 2 Consequently, (3.11) reflects the expectation that an odd derivative inμ should vanish at vanishing z.

Analytic continuation
We still cannot compute K i,j -observables defined above directly on the lattice due to the imaginary part of the action. This can be cured by analytic continuation of z. It was found in [1] that this way of analytic continuation is highly superior to a more straight-forward continuation in the full quark chemical potential µ, since the leading contribution of µ is quadratically through y. An analytic continuation in z has the advantage of leaving y unchanged and only modifying the imaginary part of the action. As for the second derivatives, the condensates contributing to the higher order derivatives also only mildly vary with iz and therefore can be expanded in powers of iz. Depending on how many powers of Φ 3 − Φ 3 a condensate contains it contains either exclusively even or exclusively odd powers of iz. Consequently, we apply the following procedure for the analytic continuation of the n-th derivative: 1. We expand the condensates with the highest sum of powers of Φ 2 and Φ 3 to first non-trivial order.
2. The coefficients of condensates with the highest power sum also re-occur in condensates with a lower power sum. We expand the lower condensates to consistent order in iz, i.e. the order where the highest power sum coefficients become relevant.

4.
We re-phrase the expanded expressions in terms of the respective condensates at imaginary chemical potential z I . As a matter of choice, we set z R = z I in order to make the factors of z I /z R unity.
As a first instructive example, we start with the analytic continuation of (3.8), as already done in [1]. We start with expanding where coeffcients a ij,k are assumed to be independent of z. In the next step, we notice that K 2,2 and K 1,1 are related by , (3.13) so in order to obtain a nontrivial contribution by K 2,2 which also contributes everywhere else consistently, we expand K 1,1 = Φ 2 ≈ a 11,0 + a 11,2 z 2 (3.14) and identify When performing the analytic continuation, we therefore end up with where all condensates K now have z I as an argument as opposed to z R before. As a matter of choice, we decide to set z R = z I to force all ratios of z R /z I to unity.
Applying the same procedure to the third derivative in (3.11), we obtain the condensates After analytic continuation we obtain The procedure outlined in this section can be easily generalized to even higher orders. We present the analytically continued expressions for the quartic, quintic, and sextic derivatives in Appendix B for the sake of readability.

Lattice implementation
Discretizing (2.1) on a spatial grid with finite lattice spacing a yields with the spatial link variables U i (x) connecting lattice sites x and x + aî, and the rescaled lattice version of the scalar field Φ L . Both x and y receive multiplicative and additive renormalization. Since EQCD is a super-renormalizable field theory, only a finite number of divergent counterterms suffice to render all amplitudes of the theory finite. In particular, only the scalar mass term y receives a divergent contribution, which can be eliminated analytically by a calculation of δy in lattice perturbation theory [31]. Beyond the elimination of divergent bits, it is also possible to pursue the lattice improvement procedure further in order to bring our lattice implementation of EQCD even closer to the continuum. While the elimination of O(a) errors in the EFT parameters of EQCD was recently completed [32] for z = 0, the O(a) behavior of most of the condensates that we measure is not yet investigated. Moreover, the emergence of the z-term in the EQCD action can also modify the known O(a) behavior of the EFT parameters. Therefore, we use a partly-O(a)-improved set of the parameters [31], meaning that β, Z 4 , δx, and Z 2 are accurate to O(a), while we set Z 3 = 1 and δy is only utilized to O(ln a), meaning that both z and y still give rise to O(a)-correction on the Lagrangian level. A compendium of the lattice renormalization constants can be found in Appendix A of [33]. Note that z does not acquire additive renormalization due to the Φ 3 → −Φ 3 -symmetry of the Lagrangian. One can argue on dimensional grounds that only K 1,1 and K 2,3 can receive divergent lattice corrections. In fact, they have been calculated [26,34] and read c 2 ≈ 0.6796(1) + ln 6 (4.7) All other condensates follow Given the very mild dependence of the condensates, especially the higher cumulants, on the lattice cutoff a, full removal of all O(a)-lattice effects is not crucial for our work (see Fig. 1) -especially of the condensates for the higher cumulants -on the lattice cutoff a, a The line in the (x, y)-plane on which EQCD describes real-world QCD lies in the supercooled phase of EQCD [26]. However, the metastability of EQCD diminishes as one approaches the crossover temperature in full QCD. Since our smallest temperature in Table 1 is too close to the crossover temperature, we applied a bit of mass reweighting to the z = 0.0, 0.025, 0.05 points at T = 277 MeV to ensure that the simulation remains in the symmetric phase. This procedure additionally inflates the errors of these points.
Our lattice EQCD simulation program is modified from the one used in [32] and is based on the openQCD-1.6-package [35] by Martin Lüscher. A combined update of 2 heatbath sweeps followed by 9 over-relaxation sweeps in a checkerboard ordering through the entire volume is applied. The scalar kinetic term in (4.1) is incorporated into the gauge sector update via a Metropolis step. The quartic scalar interaction also requires an accept/reject step in the scalar heatbath update. Both heatbath and over-relaxation are concluded by a final Metropolis step for the cubic term.
In order to extrapolate our results to the continuum, we simulate each of the six temperatures in full QCD from Table 1 times at six different values for z and at five different lattice spacings each: In total this makes 180 distinct simulations. The simulation parameters and the respective reached statistical power are tabulated in Appendix A. Due to the full O(a)-improvement in EQCD not being known in our setup, we combine the condensates to the different susceptibilities at finite a and extrapolate the susceptibilities to the continuum instead of the single condensates. As long as there are more condensates involved into the computation of a susceptibility than different variants of the susceptibility exist, this procedure is preferred since it minimizes the number of interpolations. We found it to be sufficient to fit a quadratic curve for each susceptibility, so that each fit has five data points and three fit parameters. An example for the continuum extrapolation is displayed in Fig. 1. We see from here that the overall dependence on the lattice cutoff is mild and the data can be well described with a quadratic curve in g 2 3d a. Since EQCD possesses a mass gap, finite volume effects are exponentially suppressed. Therefore, it is sufficient to choose a volume that satisfies the empirical bound [36] V a 3 ≥ To this end, we keep the physical volume fixed at (4.13)

Results
We provide our continuum-extrapolated, tabulated results in three-dimensional units in Appendix C. Comparing to four-dimensional predictions from the lattice and perturbation theory requires an application of the matching procedure outlined in Sec. 2.

Second derivatives
We start reviewing our simulation results with the qualitatively already known data for the quadratic quark number susceptibility, evaluated at different n f than in [1]. The quadratic susceptibilities are plotted as a function of the temperature in Fig. 2. For the diagonal susceptibility, we essentially observe an exponential dependence on the temperature which is generated by the exponential temperature dependence of the screening mass term y through the scale setting. At lower temperatures close to T c , the diagonal susceptibility becomes slightly negative, a feature that is already known from [1]. Increasing the chemical potentialμ results in a mild increase of χ ii 3d , the curves are strictly ordered by their respective z. As expected, the off-diagonal quark number susceptibility χ ij 3d is numerically suppressed by about one order of magnitude compared to its diagonal counterpart χ ii 3d . The off-diagonal susceptibility vanishes atμ = 0 for temperatures far away from the pseudocritical temperature of QCD. Coming closer to the crossover, the quality of the perturbative dimensional reduction decreases, so a deviation from the expected value of 0 can be explained as an artifact of a poor dimensional reduction. Deviations from 0 set in as one turns on the chemical potential, also preserving the hierarchy in z observed in the diagonal susceptibility. An important benchmark for our data is the comparison to full QCD lattice data [37] on the small-T side and full QCD perturbation theory on the large-T side in Fig. 3. To this end, we have to apply the matching procedure outlined on Sec. 2. Since we use a different number of massless quark flavors n f for different temperatures, it does not seem sensible to draw the perturbative prediction as a line. Therefore, we computed the values at the same temperatures as our lattice EQCD data. Comprehensive lattice data for all conserved charges (B, I, Q, and S) is only available for the quadratic susceptibilities, and we analyze their behavior before moving on to the higher cumulants. The conserved-chargesusceptibilities are related to our diagonal and off-diagonal quark number susceptibilities via We see that agreement with lattice data is not even reached at T = 400 MeV, although the gap between our values and the full-QCD lattice prediction shrinks by a good margin as temperature increases. This seems to contrast the earlier experience from [36]. However, there are two important differences between the present work and [1]. On the one hand, we use a MS-scale of Λ MS = 341 MeV from 2 + 1-flavor lattice QCD in contrast to the much lower Λ MS = 245 MeV from 2-flavor lattice QCD in [1]. Therefore, our smallest temperature features a much lower T /Λ MS -ratio than the corresponding temperature in [1] and consequently lies much closer to the QCD crossover, where dimensional reduction ceases to work. Moreover, our computation includes three massless flavors of quarks instead of two. At small temperatures, neglecting the strange quark mass might actually not be a valid approximation and causes a power series of extra corrections in m s /T . As expected, our results approach the perturbative result with increasing temperature and good agreement with the perturbative result is reached at T = 25 GeV, which does not contradict earlier results [1]. As the dimensional reduction technique itself is expected to work beyond T 2T c [38], the conclusion must be that nonperturbative physics plays an important role in EQCD for temperatures below T = 25 GeV. Another peculiar feature of fig. 3 is that the purely perturbative predictions are closer to the fourdimensional lattice results than our lattice EQCD results matched to full QCD. We suspect that this is due to the using a 1 loop running coupling for the matching to EQCD in (2.3), whereas [14] features a 2 loop running coupling. While there should not be any difference formally since the matching at O(g 4 ) is incomplete, anyways, there is a big difference in the numerical values of g 2 at low temperatures; the 2 loop g 2 is about 20% smaller than its 1 loop equivalent at T = 277 MeV.
Both features can be found in Fig. 4, consistency with zero for the z = 0 data is indeed realized at a quite exact level. Moreover, the derivatives in three-dimensional units contain an overall factor of π fromμ = µ πT , while the more natural scaling would be µ T . Therefore, we expect the overall order of magnitude of the third derivative to be higher than the second one, due to this extra factor of π, also translating to the ratios of the other higher cumulants in three-dimensional units. Fig. 4 furthermore shows clearly that the error of data points closer to T c drastically increases, a feature for which the small mass reweighting of some of these data points cannot be entirely blamed.
While it becomes evident from the data that there is a certain lower bound in temperature below which the dimensional reduction is not valid any more, the situation for large z is somewhat different. Even though a value of z = 0.1 at n f = 3 corresponds to µ ≈ T , which is definitely outside of the range of validity of our dimensional reduction procedure, there is no inherent evidence from the data for the breakdown of the µ Tapproximation that went into the scale setting in (2.3). Where temperatures close to the pseudocritcal temperature of full QCD cause the error bars to blow up, a large value of z just strongly favors one of the two possible signs of the Tr Φ 3 operator in the EQCD action (2.1) and therefore reduces the fluctuations in the system, leading to smaller error bars at comparable numerical effort.

Fourth derivatives
The fourth derivative is again even inμ, therefore, the fully diagonal susceptibility χ iiii 3d and its bi-diagonal counterpart χ iijj 3d do have a finite expectation value even atμ = 0, which also reflects in our results in Fig. 5. Despite the lack of a clear physical interpretation, we still display χ ijkl 3d at temperatures with n f > 3, since it gives an idea about the part of (B.1) that solely contains condensates that generically appear in fourth derivatives. A rough scaling by an additional factor of π compared to the third derivatives is maintained here, too. Following the trend of the third derivatives, the data at T = 277 MeV starts to become ever more noisy, whereas the error bars for all other temperatures are barely visible.

Fifth derivatives
Fifth derivatives are again odd inμ and therefore expected to feature a vanishing expectation value atμ = 0, as one can see from our results in Fig. 6. Observing another relative increase of π in the order of magnitude, the consistency with 0 of the z = 0-data is still very precise. Aside from the generically fifth-derivative-contribution χ ijklm

3d
, Fig. 6 only shows data from scenarios in which the respective derivative is well-defined, i.e. χ iijkl 3d is only displayed at temperatures that correspond to at least n f = 4. Even though without physical meaning, we still display the numerical values of the derivatives omitted in the plot in App. C.

Sixth derivatives
Just as for the previous derivatives, we observe a relative factor of π in the order of magnitude of the sixth derivative. The relative errors of T = 277 MeV increased once more, only physically well-defined derivatives are plotted other than the fully off-diagonal part χ ijklmn

Comparison of higher derivatives to perturbation theory
Beyond second derivatives, comparing our data for the higher derivatives to their respective perturbative predictions provides an important benchmark. Restricting ourselves tō µ = 0, we plot all susceptibilities which are not expected to vanish for symmetry reasons in Fig. 8. Just like for the second derivatives, the discrepancy near the QCD crossover temperature is substantial. As one increases the temperature, the two curves approach each other. The lowest temperature was not even included in the plots of the sixth order derivatives since it essentially consists of a huge error bar and would only prevent the reader from noticing details in the plot at higher temperatures. The higher derivatives feature a quicker convergence towards the perturbative solution in general. We translate our quartic (e) (f) Figure 8. Comparison of perturbative results [14], hotQCD lattice results [21], and our results for quartic and sextic susceptibilities atμ = 0.
susceptibilities to derivatives with respect to the baryon chemical potential µ B via Technically, the above expression is only exact for n f = 3. However, if applying it consistently to perturbative results and EQCD results, both relying on the same number of massless flavors, our data and perturbation theory are still comparable, even at higher n f . Fig. 8(f) compares our data to perturbation theory and lattice data [21]. Since the lattice data is only available at two different N t , an extrapolation to the continuum is not possible. However, we face the same significant deviations of our low-temperature results from lattice computations. The margin is so big here that it is excluded that this can be tracked down to remaining lattice-cutoff effects of the full QCD data. On the contrary, the reasons for this deviation seem to be the same as for the big deviations from the lattice data for the second derivative: neglecting the strange quark mass, poor scale fixing, and insufficient accuracy in the matching of the effective theory close to the pseudocritical temperature.

Conclusion
EQCD is an effective description of high-temperature QCD which can be robustly derived using perturbation theory. It fully includes the non-perturbative soft thermal physics, and offers an economical way to study its effects in Monte Carlo simulations which are numerically less costly than full 4-dimensional lattice simulations. It also has a dramatically milder sign problem than full QCD. In the present work, we have derived formulae for higher derivatives of the pressure with respect to the quark chemical potential up to sixth order. We have studied these operators on the lattice and provided continuum-extrapolated results in six cases corresponding to physical scenarios at different temperatures T and number of massless quark flavors n f . Our data was compared to perturbation theory at high temperatures and full-QCD lattice simulations at low temperatures close to the crossover temperature T c of full QCD at vanishing quark chemical potential. For temperatures reasonably well above T c , we were able to extract a clear signal on the lattice for all derivatives up to sixth order; a clean extrapolation to the continuum is possible. Agreement with perturbative results sets in around T ≈ 25 GeV.
Our results shed light into the limitations of the dimensional reduction procedure in general, and the perturbative treatment of EQCD. The data at n f = 3 show substantial deviations from n f = 2 + 1 full QCD lattice data. We conclude that the dimensional reduction procedure at O(g 3 )-accuracy does not work for temperatures below 400 MeV, probably even higher temperatures. In particular, the reasons for this are threefold: First, the quality of the dimensional reduction is poor around the pseudocritical temperature of full QCD, which is where full QCD-lattice data is available. Especially the matching of the chemical potential-term of full QCD to the three-dimensional z constrains the accuracy of the matching to O(g 3 ). The substantially better agreement of purely perturbative results (containing EQCD treatment for the soft sector) with fourdimensional lattice predictions hints that already a -formally incomplete -usage of a higher loop running coupling could accelerate the convergence of our fully matched lattice EQCD results towards fourdimensional lattice predictions drastically. Second, we neglected the mass of the strange quark, which turns out not to be a valid approximation up to T ≈ 400 MeV. A finite strange quark mass could be incorporated into the dimensional reduction procedure from full QCD to EQCD [39], even though that would partly break the symmetry of the derivatives in flavor space. Third, and somewhat related to the previous point, we set the scale by a value for Λ MS for 2 + 1 flavor full QCD, which assumes a finite strange quark mass and therefore deviates from the n f = 3-flavor results that would be more rigorous to use in our case.
Under the caveat of these improvements, lattice EQCD simulations of higher order cumulants could deliver a valuable tool for the determination of the equation of state of nuclear matter at temperatures of a few times the pseudocritical temperature and moderate quark chemical potential.

Acknowledgments
We thank Aleksi Vuorinen and Mikko Laine for fruitful discussions on the present work, in particular the perturbative aspects. We highly appreciated Olaf Kaczmarek and Frithjof Karsch providing us with raw data for the higher cumulants from the lattice. Moreover, we thank Ari Hietanen for his initial work on the analytical expressions of the third and fourth derivatives from EQCD. The support of the Academy of Finland grants 308791 and 320123 is acknowledged. Table 2 shows the simulation parameters at which we conducted out three-dimensional Monte-Carlo simulations. A total computation budget of 430, 000 CPU × hrs was utilized. Due to bad core allocation on the local cluster, some of the statistics of different temperature or z, but equal g 2 3d a, differ a bit. As already pointed out, we applied a small amount of mass reweighting to the z = 0.0, 0.025, 0.05 data points at T = 277 MeV at all considered lattice spacings to stabilize the system in the supercooled phase. The reweighting allowed a computation at y comp = 0.375.

B Expressions for higher derivatives
In this appendix, we present the full analytically continued expressions for the higher derivatives. Each new derivative requires new condensates and new way of analytically continuing the lower condensates. We denote the condensates with K. The first index is the sum of the two powers of Φ 2 and Φ 3 . The second index runs from 1 to firstindex + 1 and refers to the power of Φ 3 involved. Symmetrization in indices are denoted by three dots.

B.1 Fourth derivatives
Continuing the condensates analytically yields the analytically continued expression for the fourth derivative We see that K 3,2 is the first derivative of K 2,1 with respect to iz, K 4,2 the first derivative of K 3,1 , and K 4,4 the first derivative of K 3,3 . Therefore, the three named condensates of lower order behave non-trivially under the analytic continuation.

B.2 Fifth derivatives
Analytically continuing the condensates as

B.3 Sixth derivatives
The sixth order derivative contains the condensates which can be combined to the analytically continued expression for the sixth order deriva- Beyond the usual identification of condensates that correspond to first or second derivatives of other condensates with respect to iz, we see the connection of three condensates for the first time via

C Tabulated results
In this appendix, we give tabulated numerical results for all (three-dimensional) cumulants that we computed, in units of the three-dimensional coupling g 2 3d . The number of different indices of some of the derivatives exceeds the number of massless quark flavor available. Therefore, these cumulants do not correspond to a physical scenario. However, χ ijkl 3d for instance has no physical meaning in scenarios with n f = 3, but gives a good estimate how big the influence of the generically fourth-order condensates in (B.1) is, which still contains valid, though not strictly physical, information.    Table 5. Results for the different derivatives of the EQCD pressure at T = 25, 100 GeV, number of massless quark flavors n f , and chemical potentials z = n f µ 3π 2 T ; in units of the three-dimensional coupling g 2 3d .