Three-loop HTL QCD thermodynamics

The hard-thermal-loop perturbation theory (HTLpt) framework is used to calculate the thermodynamic functions of a quark-gluon plasma to three-loop order. This is the highest order accessible by finite temperature perturbation theory applied to a non-Abelian gauge theory before the high-temperature infrared catastrophe. All ultraviolet divergences are eliminated by renormalization of the vacuum, the HTL mass parameters, and the strong coupling constant. After choosing a prescription for the mass parameters, the three-loop results for the pressure and trace anomaly are found to be in very good agreement with recent lattice data down to $T \sim 2-3\,T_c$, which are temperatures accessible by current and forthcoming heavy-ion collision experiments.


Introduction
The last generation of relativistic heavy-ion collision experiments exceeded the energy density necessary for the formation of a quark-gluon plasma, motivating the development of a quantitative framework with which to calculate the properties of this new state of matter. At Brookhaven National Laboratory (RHIC), New York, USA, the initial temperatures were up to twice the critical temperature for deconfinement, 1 , T c ∼ 170 MeV. This translates to a strong coupling constant of α s (µ = 2π × 170 MeV) ≡ g 2 /(4π) ∼ 0.43 where µ is the renormalization scale and it relates to the temperature by µ = 2πT . The upcoming experiments at CERN (LHC), Geneva, Switzerland, are expected to yield initial temperatures of 4−6 T c , driving the running coupling further down. Initially, the hope was that the asymptotic freedom of QCD would allow calculations using a perturbative expansion in the coupling. Utilizing a weak-coupling expansion in the coupling constant g to calculate the thermodynamic functions of QCD has a long history [2][3][4][5][6][7][8][9], and the pressure is now known through order g 6 log g for non-Abelian gauge theories [10]. Unfortunately, for all but tiny coupling constants, and thus astronomically high temperatures, these expansions converge very poorly, and show large dependence on the renormalization scale. In figure 1, we show the weak-coupling expansion for the QCD pressure with N f = 3 normalized to that of an ideal gas through order g 5 . The various approximations oscillate wildly and show no signs of convergence in the temperature range shown. The bands are obtained by varying the renormalization scale µ by a factor of 2 around the value µ = 2πT and we use three-loop running for α s [11] with Λ MS (N f = 3) = 344 MeV [12]. This oscillating behavior is generic for hot field theories, and not specific to QCD. The instability is thought to be caused by plasma effects such as screening of electric fields and Landau damping. This calls for a nonperturbative approach, or a reorganization of the perturbative expansion that takes such effects into account. Furthermore, data from RHIC suggest that the matter created behaves more like a strongly coupled fluid with small viscosity [13][14][15][16][17], inspiring the development of strongly coupled formalisms, perhaps the most successful being those based on the Anti-de-Sitter space/conformal field theory (AdS/CFT) correspondence proposed by Maldacena [18]. However, work on the perturbative side has shown that observables like jet quenching [19,20] and elliptic flow [21] can also be described by a perturbative setup. Looking forward to the upcoming heavy-ion experiments scheduled to take place at the LHC, it is therefore important to know if, at the higher temperatures generated, one should expect a strongly-coupled (liquid) or weakly-coupled (plasma) description to be more appropriate. The key question is, will the generated matter behave more like a plasma of quasiparticles at these higher temperatures.
Of course, another approach is lattice gauge theory, which is a nonperturbative discrete space-time framework that is the closest one can currently get to a first-principles calculation with realistic parameters. However, the Monte Carlo methods used currently restrict lattice gauge theory to the study of systems with low chemical potential, and reliable calculations of dynamical quantities have proved difficult. A quark-gluon plasma created in a heavy-ion collision will be out of thermal equilibrium and will have nonzero baryon density, making these restrictions quite severe. Still, the part of the phase diagram that is accessible by lattice gauge theory becomes a clean testing ground for the various approaches to analytical QCD calculations. In this paper we will be comparing our new "reorganized" perturbative results for the thermodynamic functions to recent lattice data from refs. [22] and [23].
There are several ways of systematically reorganizing the perturbative expansion [24][25][26]. In screened perturbation theory (SPT) [27][28][29][30][31] which was inspired by variational perturbation theory (VPT) [32][33][34][35][36][37], a mass term is added to the free Lagrangian, and then subtracted as an interaction term. Unfortunately, this reorganization cannot be directly applied to gauge theories, since the added local mass term would violate gauge invariance [38][39][40]. Instead, one can add and subtract a hard-thermal-loop (HTL) effective action [41] which modifies the propagators and vertices systematically and self-consistently, so that the reorganization is manifestly gauge invariant [42][43][44]. The mass parameters m D and m q are introduced, and identified with the Debye screening mass and the thermal quark mass at leading order, respectively, in order to reproduce the high temperature limit of thermal QCD; however, at higher orders it is necessary to introduce a prescription for the subleading contributions to the mass parameters. The resulting HTL perturbation theory (HTLpt) framework can be applied to static as well as dynamic quantities.
The HTLpt framework has recently been applied to QED [45], where the calculation of the thermodynamic functions is carried out through next-to-next-to-leading-order (NNLO). The thermodynamic functions of QCD have been calculated to next-to-leading-order (NLO) [46,47]. For pure-glue QCD the thermodynamic functions were recently calculated to NNLO in [48,49]. Our paper builds upon the results of NNLO QED and pure-glue QCD, and we now include NNLO contributions from quark and quark-gluon interaction diagrams to express the thermodynamic functions of full QCD to NNLO. The calculation is presented in some detail. For a short letter with just the final results, we refer to ref. [50]. Our results indicate that the lattice data at temperatures T ∼ 2 T c are consistent with the quasiparticle picture. This is a nontrivial result since, in this temperature regime, the QCD coupling constant is neither infinitesimally weak nor infinitely strong with g ∼ 2, or equivalently α s ∼ 0.3. Therefore, we have a crucial test of the quasiparticle picture in the intermediate coupling regime.
After the introduction, we begin in section 2 with a summary of HTLpt applied to QCD. Section 3 contains expressions for the diagrams contributing to the thermodynamic potential up to NNLO. In section 4 these diagrams are reduced to scalar sum-integrals and expanded in powers of m D /T and m q /T , keeping all terms that contribute through order g 5 if the mass parameters are taken to be of order g at leading order (LO). In section 5 the calculated diagrams are gathered to obtain the renormalized NNLO thermodynamic potential. The results are then presented in section 6, and compared to recent lattice data. Due to difficulties with the normal variational approach, there is special emphasis placed on the determination of the mass parameters m D and m q . In section 7 we summarize.

HTL perturbation theory
The Lagrangian density for QCD in Minkowski space is where the gluon field strength is , the term with the quark fields ψ contains an implicit sum over the N f quark flavors, and the covariant derivative is D µ = ∂ µ − igA µ . The ghost term L gh depends on the gauge-fixing term L gf . In this paper we choose the class of covariant gauges where the gauge-fixing term is The perturbative expansion in powers of g generates ultraviolet divergences. The renormalizability of perturbative QCD guarantees that all divergences in physical quantities can be removed by renormalization of the coupling constant α s = g 2 /(4π) and the necessary counterterms are represented by ∆L QCD in the Lagrangian (2.1). There is no need for wavefunction renormalization, because physical quantities are independent of the normalization of the field. There is also no need for renormalization of the gauge parameter, because physical quantities are independent of the gauge parameter.
Hard-thermal-loop perturbation theory (HTLpt) is a reorganization of the perturbation series for thermal QCD. The Lagrangian density is written as The HTL improvement term is where y µ = (1,ŷ) is a light-like four-vector, and . . . y represents the average over the directions ofŷ. The term (2.4) has the form of the effective Lagrangian that would be induced by a rotationally invariant ensemble of charged sources with infinitely high momentum and modifies the propagators and vertices self-consistently so that the reorganization is manifestly gauge invariant [38,[51][52][53][54][55]. The parameter m D can be identified with the Debye screening mass, and m q with the thermal quark mass to account for the screening effects. HTLpt is defined by treating δ as a formal expansion parameter. By coupling the HTL improvement term (2.4) to the QCD Lagrangian (2.1), HTLpt systematically shifts the perturbative expansion from being around an ideal gas of massless particles which is the physical picture of the weak-coupling expansion, to being around a gas of massive quasiparticles which are the more appropriate physical degrees of freedom at high temperature. Physical observables are calculated in HTLpt by expanding them in powers of δ, truncating at some specified order, and then setting δ = 1. This defines a reorganization of the perturbation series in which the effects of m 2 D and m 2 q terms in (2.4) are included to all orders but then systematically subtracted out at higher orders in perturbation theory by the δm 2 D and δm 2 q terms in (2.4). If we set δ = 1, the HTLpt Lagrangian (2.3) reduces to the QCD Lagrangian (2.1).
If the expansion in δ could be calculated to all orders the final result would not depend on m D and m q when we set δ = 1. However, any truncation of the expansion in δ produces results that depend on m D and m q . Some prescription is required to determine m D and m q as a function of T and α s . We will discuss several prescriptions in section 6.
The HTL perturbation expansion generates ultraviolet divergences. In QCD perturbation theory, renormalizability constrains the ultraviolet divergences to have a form that can be cancelled by the counterterm Lagrangian ∆L QCD . We will demonstrate that the renormalization of HTLpt can be implemented by including a counterterm Lagrangian ∆L HTL among the interaction terms in (2.3). There is no proof that the HTL perturbation expansion is renormalizable, so the general structure of the ultraviolet divergences is not known; however, it was shown in previous papers [46,47] that it was possible to renormalize the next-to-leading order HTLpt prediction for the pressure of QCD using only a vacuum counterterm, a Debye mass counterterm, and a fermion mass counterterm. By also including a coupling constant counterterm, this is possible at NNLO as well, as we will show in this paper. The necessary counterterms for the renormalization at NNLO as just discussed read where the vacuum and thermal mass counterterms were derived in ref. [46,47], and the coupling constant counterterm is the standard one-loop running for QCD [56,57].
One-loop pure gauge contribution to gluon self-energy Figure 2. QCD diagrams contributing through NLO in HTLpt. The spiral lines are gluon propagators, the dotted lines are ghost propagators, and the solid lines are quark propagators. A circle with a Π indicates a one-loop gluon self-energy insertion, a Σ indicates a one-loop quark self-energy insertion, and a Γ indicates a one-loop vertex insertion. A square with a g is shorthand for the pure gauge diagrams contributing to the one-loop gluon self-energy. All gluon and quark propagators and vertices shown are HTL-resummed propagators and vertices. The logic behind the diagram notation is as follows: diagrams consisting only of gauge propagators have g superscripts. Diagrams containing fermion propagators have f superscripts. The subscript indices are identical to those used in [49] (pure gauge QCD) and [45] (QED). We do not display the symmetry factors in the diagrams. Figure 3. QCD diagrams contributing to NNLO in HTLpt.

Diagrams for the thermodynamic potential
Although HTLpt is defined in Minkowski space, it is much more convenient to carry out the necessary calculations in Euclidean space since in the paper we focus on static quantities. In the imaginary-time formalism, Minkoswski energies have discrete imaginary values p 0 = i2nπT for bosons and p 0 = i(2n + 1)2πT for fermions, and integrals over Minkowski space are replaced by sum-integrals over Euclidean vectors (2nπT, p) or (2(n + 1)πT, p), respectively, with n an integer. We will use the notation P = (P 0 , p) for Euclidean momenta. The magnitude of the spatial momentum will be denoted p = |p|. The inner product of two Euclidean vectors is P · Q = P 0 Q 0 + p · q. The vector that specifies the thermal rest frame remains n = (1, 0). Here we use dimensional regularization and the remaining three-dimensional integral is generalized to d = 3 − 2 spatial dimensions. We define the dimensionally regularized sum-integral by where µ is an arbitrary momentum scale. The factor (e γ E /4π) is introduced so that, after minimal subtraction of the poles in due to ultraviolet divergences, µ coincides with the renormalization scale of the MS renormalization scheme. The Feynman diagrams discussed in this section are gathered in figures 2 (notation key, LO, NLO) and 3 (NNLO). Using the same notation for the group theory factors as Arnold and Zhai [6], for QCD with N c colors and N f flavors of quarks, we have 3) The earlier work on NNLO pure-glue QCD [49] and NNLO QED [45] essentially treat special cases where the group symmetry factors are and as such the analytical results from those papers will be highly useful in the present calculation as well.
The thermodynamic potential at LO in HTLpt for QCD reads Here, F g 1a is the contribution from the gluons and F g 1b the contribution from the ghost, and ∆ 0 E 0 is the leading-order vacuum counterterm. The expression of eq. (3.6) is presented in ref. [47], and the details can be found in section 3 therein.
The thermodynamic potential at NLO reads where ∆ 1 E 0 , ∆ 1 m 2 D , and ∆ 1 m 2 q are the terms of order δ in the vacuum energy density and mass counterterms. Again the detailed expression of eq. (3.7) is presented in section 3 of ref. [47].

Expansion in the mass parameters
In refs. [46,47], the NLO HTLpt pressure was reduced to scalar sum-integrals. It was clear that evaluating these scalar sum-integrals exactly was intractable and the sum-integrals were calculated approximately by expanding them in powers of m D /T and m q /T following the method developed in ref. [30]. We will adopt the same strategy in this paper and carry out the expansion to high enough order to include all terms through order g 5 if m D and m q are taken to be of order g. The pressure can be divided into contributions from hard and soft momenta, which are the momenta proportional to the scales T and gT respectively. In the one-loop diagrams, the contributions are either hard (h) or soft (s), while at the two-loop level, there are hard-hard (hh), hard-soft (hs), and soft-soft (ss) contributions. At three loops there are hard-hard-hard (hhh), hard-hard-soft (hhs), hard-soft-soft (hss), and soft-soft-soft (sss) contributions. As mentioned above, the (h), (s), (hh), (hs) and (ss) contributions have been calculated in refs. [46,47], while most of the three-loop contributions can be read out from the earlier NNLO HTLpt work on pure-glue QCD [49] and QED [45], thus to keep the discussion compact, we only list the contribution that is specific to QCD here, i.e the term proportional to c A s F in eq. (3.8).
If all momenta in the three loops are hard, we can obtain the m D /T and m q /T expansion by expanding in powers of m 2 D and m 2 q , and in this case to obtain the expansion through order g 5 , we can simply use bare propagators and vertices. The contributions from the three-loop diagrams were first calculated by Arnold and Zhai in ref. [5], and later by Braaten and Nieto in ref. [8] from which the (hhh) contribution to the c A s F term reads Note that there is no (hhh) contribution from the HTL-specific diagram F f 3o [45]. In the (hhs) region, the momentum P is soft, 2 while the momenta Q and R are always hard. The function that multiplies the static propagator ∆ T (0, p), ∆ L (0, p) or ∆ X (0, p) can be expanded in powers of the momentum p. In the case of ∆ T (0, p), the resulting integrals over p have no scale and thus vanish in dimensional regularization. The integration measure p scales like m 3 D , the static propagators ∆ L (0, p) and ∆ X (0, p) scale like 1/m 2 D , and every power of p in the numerator scales like m D . 3 The (hhs) contribution to the c A s F term reads The soft contribution arises from the P0 = 0 term in the sum-integral. See appendix A of ref. [58] for details. 3 We refer the reader to e.g. ref. [46] for HTLpt Feynman rules, notations and conventions.
with the function T P defined by where the angular brackets represent an average over c For all of the diagrams that are infrared safe, the (hss) contribution is of order g 4 m 2 D , i.e. g 6 , and can be ignored. The infrared divergent diagrams, i.e. F f 3n , contribute to the c A s F term as follow (4.5)

The thermodynamic potential
In this section we present the final renormalized QCD thermodynamic potential. The LO and NLO results were derived in refs. [46,47], thus to keep the presentation compact, we only list the NNLO one, which corresponds to order δ 2 in HTLpt.
Using the results listed in section 4, the renormalization contributions at order δ 2 read where F ideal is the pressure of a gas of d A massless spin-one bosons andm D ,m q andμ are dimensionless variables: In order to keep the expression (5.1) compact, the pure-glue contributions have been represented by ∆Ω YM 2 which is given by eq. (4.10) of ref. [49]. Adding the NNLO counterterms (5.1) to the contributions from the various NNLO diagrams we obtain the renormalized NNLO thermodynamic potential. We note that at NNLO all numerically determined subleading coefficients in drop out and due to this the final result is completely analytic. The resulting NNLO thermodynamic potential is where Ω YM NNLO is the NNLO pure-glue thermodynamic potential given by eq. (4.11) in ref. [49] We note that the coupling constant counterterm listed in eq. (2.8) coincides with the known one-loop running of QCD In the next section we will present results as a function of g evaluated at the renormalization scale 2πT .

Thermodynamic functions
The mass parameters m D and m q in HTLpt are completely arbitrary. To complete a calculation, it is necessary to specify m D and m q as functions of g and T . As HTLpt is inspired by variational perturbation theory, the natural prescription for specifying the mass parameters is to use the variational solutions resulting from the gap equations The LO gap equations are not well defined due to the lack of the coupling constant g in the LO thermodynamic potential, and in previous works on LO HTLpt one therefore chose the LO weak-coupling perturbative values of the Debye and thermal quark masses as physical estimates of the HTLpt mass parameters [42][43][44]. At NLO the variational prescription gives nontrivial solutions for the mass parameters, and the NLO results presented in [46,47] use this prescription for both m D and m q . For our NNLO results the variational prescription yields a complex Debye mass and m q = 0. The problem of a complex mass parameter was encountered at NNLO in SPT massless φ 4 theory [29], in HTLpt QED [45], and in HTLpt pure-glue QCD [48,49], so this seems to be a general issue with SPT/HTLpt at NNLO, which is not specific to QCD. Since the weak-coupling perturbative Debye mass receives contributions from the nonperturbative magnetic scale beyond LO, for the NNLO HTLpt pure-glue QCD results the Braaten and Nieto's (BN) mass parameter of three-dimensional electrostatic QCD (EQCD) was used as a substitute for the perturbative Debye mass, effectively discarding the nonperturbative contributions. The perturbative quark mass does not suffer from this issue, and as such there is bound to be some "mixing" of prescriptions for the mass parameters of the NNLO full QCD results unless one accepts the variational complex m D value and sticks to the variational prescription for both parameters. In this section we will discuss and compare numerically the aforementioned prescriptions, before making a choice of prescriptions to use for comparison of our results to lattice data.

Variational masses
Solving the NNLO gap equations (6.1)-(6.2) yields a complex Debye mass and m q = 0. Figure 4 shows the real and imaginary parts of the complex NNLO Debye mass, scaled by the LO perturbative m D . While the imaginary contribution is non-negligible at intermediate coupling, and even grows to surpass the real contribution as α s approaches unity, it disappears in the small coupling limit, where the NNLO variational and LO perturbative results for m D overlap. Presumably the complex mass parameter is an artifact of the truncation after the m D /T and m q /T expansions, but this is difficult to confirm without extending the truncation to higher order. One strategy is then to throw away the imaginary part of the thermodynamic potential to obtain thermodynamic functions that are real valued. Figure 5 shows the NNLO pressure found using this approach. For reference, figure 6 compares the real valued and the discarded imaginary contributions to the NNLO variational pressure over the same  temperature range.

Perturbative masses
At LO in the coupling constant g, the Debye mass is given by the static longitudinal gluon self-energy at zero three-momentum, m 2 D = Π L (0, 0), i.e.
At NLO the weak-coupling perturbative Debye mass becomes logarithmically infrared divergent, reflecting the contribution from the nonperturbative magnetic scale g 2 T [59, 60]. The nonperturbative contribution was calculated in ref. [61] and reads where m mag is the nonperturbative magnetic mass. Using this mass prescription for the Debye mass would then require input from e.g. lattice simulations. For this reason we will not consider this prescription and instead will use the Braaten-Nieto prescription detailed in the next subsection. The NLO quark mass is infrared safe, and was calculated for N c = 3 and N f = 2 in ref. [62], However, to our knowledge there has not been results available for general N f , and in order to make the comparison with lattice data feasible, we use instead, whenever it applies, the N c = 3 LO quark mass m q = gT √ 6 . (6.6)

BN mass
The strategy of equating m D to the BN mass has earlier been applied to NNLO HTLpt Yang-Mills theory [48,49]. Inspired by dimensional reduction, one equates the m D mass parameter with the mass parameter of three-dimensional EQCD [7], i.e. m D = m E . This mass can be interpreted as the contribution to the Debye mass from the hard scale T and is well defined and gauge invariant order-by-order in perturbation theory. However, beyond NLO it will also depend on the factorization scale that separates the hard scale and the soft scale gT . In ref. [7] it was calculated to NLO, giving (6.7)

Choice of mass prescriptions
To avoid dealing with imaginary contributions to the thermodynamic potential, we opt for the BN prescription for the Debye mass parameter, the same choice that was made for NNLO HTLpt pure-glue QCD [49]. As for m q , there is no problem with getting zero thermal quark mass from the variational prescription: even if the m 2 q -terms of Eq. (5.6) drop out, the quarks still contribute to the HTLpt thermodynamical potential through the terms proportional to d F , s F and s 2F . On the other hand, when using the BN mass parameter one might want to use the perturbative approach for the quark mass parameter, as using BN mass and perturbative thermal quark mass can be argued to be less of a "mixing" of prescriptions than BN mass and m q = 0. As it turns out, the final NNLO results are very insensitive to whether one chooses a perturbative mass prescription for m q , or uses the variational mass m q = 0. In figure 7 the NNLO pressure with LO perturbative m q is virtually indistinguishable from the m q = 0 result. The difference in dependence on the renormalization scale is also very small, though slightly in favor of the perturbative m q , see figure 8. However, convergence is improved with the choice m q = 0, where the pressure curves of figure 7 monotonically approaching the NNLO results rather than oscillating as one goes from LO to NLO to NNLO. We will therefore use m q = 0 in the following unless otherwise stated.

Comparison to lattice data
In figure 9 we show the normalized pressure for N c = 3 and N f = 2 + 1 (left panel), and N c = 3 and N f = 2 + 1 + 1 (right panel) as a function of T . The results at LO, NLO, and NNLO use the perturbative Debye mass given by eq. (6.7) as well as m q = 0. For the strong coupling constant α s , we used three-loop running [11] with Λ MS = 344 MeV which for N f = 3 gives α s (5 GeV) = 0.2034 [12]. The central line is evaluated with the renormalization scale µ = 2πT which is the value one expects from effective field theory calculations [7,63] and the band represents a variation of µ by a factor of 2 around this value.
The lattice data from the Wuppertal-Budapest collaboration use the stout action. Since their results show essentially no dependence on the lattice spacing (it is smaller than the statistical errors), they provide a continuum estimate by averaging the trace anomaly measured using their two smallest lattice spacings corresponding to N τ = 8 and N τ = 10 Figure 10. Comparison of NNLO predictions for the scaled trace anomaly with N f = 2 + 1 (left panel) and N f = 2 + 1 + 1 fermions (right panel) lattice data from Bazavov et al. [22] and Borsanyi et al. [23]. We use N c = 3, three-loop running for α s , µ = 2πT , and Λ MS = 344 MeV. Shaded band shows the result of varying the renormalization scale µ by a factor of 2 around µ = 2πT . See main text for details. [23], which were essentially on top of the N τ = 6 measurement [64]. 4 Using standard lattice techniques, the continuum-estimated pressure is computed from an integral of the trace anomaly. The lattice data from the hotQCD collaboration are their N τ = 8 results using both the asqtad and p4 actions [22]. The hotQCD results have not been continuum extrapolated and the error bars correspond to only statistical errors and do not factor in the systematic error associated with the calculation which, for the pressure, is estimated by the hotQCD collaboration to be between 5 -10%. We note that there are hotQCD results for physical light quark masses [65]; however, these are available only for temperatures below 260 MeV and the results are very close to the results shown in the figures so we do not include them here.
As can be seen from figure 9 the successive HTLpt approximations represent an improvement over the successive approximations coming from a weak-coupling expansion; however, as in the pure-glue case [48,49], the NNLO result represents a significant correction to the LO and NLO results. That being said, the NNLO HTLpt result agrees quite well with the available lattice data down to temperatures on the order or 2 T c ∼ 340 MeV for both N f = 3 (left panel) and N f = 4 (right panel). Below these temperatures the successive approximations give large corrections with the correction from NLO to NNLO reaching 100% near T c .
In figure 10, we show the NNLO approximation to the trace anomaly (interaction measure) normalized to T 4 as a function of T for N c = 3 and N f = 2 + 1 (left panel) and for N c = 3 and N f = 2 + 1 + 1 (right panel). 5 In the left panel we show data from both the Wuppertal-Budapest collaboration and the hotQCD collaboration taken from the same data sets displayed in figure 9 and described previously. In the case of the hotQCD results we note that the results for the trace anomaly using the p4 action show large lattice size affects at all temperatures shown and the asqtad results for the trace anomaly show large lattice size effects for T ∼ > 200 MeV. In the right panel we display a parameterization (solid blue curve) of the trace anomaly for N f = 4 published by the Wuppertal-Budapest collaboration [23] since the individual data points were not published. In both the left and right panels we see very good agreement with the available lattice data down to temperatures on the order of T ∼ 2 T c .

Large N f
In the limit N f 1 while keeping g 2 N f ∼ 1, only ring diagrams contribute to the pressure in perturbation theory. Since s F is proportional to N f in both QCD and QED, this indicates the equivalence of QED and QCD in the large-N f limit, with the large N f effective coupling defined as g eff ≡ g N f /2 for QCD, and g eff ≡ e N f for QED. The NNLO HTLpt results in this limit were discussed in the context of QED in ref. [68].
It is possible to solve for the O(N 0 f ) contribution to the pressure exactly [66], enabling comparison between predictions from approximations such as perturbation theory or HTLpt, and exact numerical results. In figure 11, we plot the NLO and NNLO HTLpt predictions for the large-N f pressure along with the numerical result of [66] as well as the perturbative g 4 eff , g 5 eff , and g 6 eff predictions at µ = e −γ E πT obtained in [67]. For the pressure, going to NNLO in HTLpt extends the range of agreement with exact results compared to NLO, from g eff 2 to g eff 2.8. The HTLpt large coupling behavior qualitatively matches that of the numerical results.

Summary and outlook
We have calculated the NNLO contributions to the thermodynamic functions of SU(N c ) Yang-Mills theory with N f fermions using HTLpt. Using the BN mass for m D and m q = 0 as the mass prescriptions, at NNLO we find that HTLpt predictions for the pressure and the trace anomaly are in agreement with lattice data for N c = 3 and N f ∈ {3, 4} down to T ∼ 2 T c . The failure of HTLpt to match lattice data at lower temperatures is to be expected since one is expanding around the trivial vacuum A µ = 0 and therefore neglects the approximate center symmetry Z(N c ), which becomes essential as one approaches the deconfinement transition [69][70][71][72][73][74][75]. In addition, it is also in line with expectations since below T ∼ 2 − 3 T c a simple "electric" quasiparticle approximation breaks down due to nonperturbative chromomagnetic effects [59,60]. 5 The Wuppertal-Budapest data were obtained using a physical charm mass at low temperatures; however, we use massless quarks. The difference between massive and massless quarks is expected to be significant only for T ∼ < 414 MeV corresponding to the temperature where the lowest fermionic Matsubara mode equals the charm quark mass. Figure 11. The predictions for the large-N f pressure of QCD between the numerical exact result from [66], NLO and NNLO HTLpt, and perturbative g 4 eff through g 6 eff [67] results at µ = e −γ E πT .
We find that when including quarks the agreement with lattice data is greatly improved as compared to the NNLO results of pure-glue QCD [48,49]. Fermions are perturbative in the sense that they decouple in the dimensional-reduction step of effective field theory, so we expect that including contributions from quarks gives at least as good agreement with the lattice calculations as the pure-glue case. However, the exact reason for the better agreement between the HTLpt predictions and lattice calculations when including quarks is not clear to us.
Just as for NNLO pure-glue QCD we found that the variational solution for the mass parameter m D is complex and we therefore chose instead to use the Debye mass from EQCD for our m D parameter. In the variational prescription the quark gap equation gives m q = 0, but even if the quarks do not get a thermal mass, the thermodynamical potential still gets HTLpt contributions from the inclusion of quarks in the calculation. Whether the complexity of the variational Debye mass parameter is due to the additional expansion in m D /T and m q /T is impossible to decide at this stage. The correction to the pressure going from NLO to NNLO is also rather large. It is unfortunate that the nonperturbative magnetic scale prevents going to N 3 LO without supplementing the calculation with input from three-dimensional lattice calculations, as it would be interesting to see whether the complexity of the Debye mass parameter and the slow convergence persists.
In closing, we emphasize that HTLpt provides a gauge invariant reorganization of perturbation theory for calculating static and dynamic quantities in thermal field theory. Given the good agreement with lattice data for thermodynamic quantities down to temperatures that are relevant for LHC, it would therefore be interesting and challenging to apply HTLpt to the calculation of dynamic quantities, especially transport coefficients, at these temperatures.