Three-loop HTL gluon thermodynamics at intermediate coupling

We calculate the thermodynamic functions of pure-glue QCD to three-loop order using the hard-thermal-loop perturbation theory (HTLpt) reorganization of finite temperature quantum field theory. We show that at three-loop order hard-thermal-loop perturbation theory is compatible with lattice results for the pressure, energy density, and entropy down to temperatures $T\simeq3\;T_c$. Our results suggest that HTLpt provides a systematic framework that can used to calculate static and dynamic quantities for temperatures relevant at LHC.

The goal of ultrarelativistic heavy-ion collision experiments is to generate energy densities and temperatures high enough to create a quark-gluon plasma. One of the chief theoretical questions which has emerged in this area is whether it is more appropriate to describe the state of matter generated during these collisions using weakly-coupled quantum field theory or a strong-coupling approach based on the AdS/CFT correspondence. Early data from the Relativistic Heavy Ion Collider (RHIC) at Brookhaven National Labs indicated that the state of matter created there behaved more like a fluid than a plasma and that this "quark-gluon fluid" is strongly coupled [1].
In the intervening years, however, work on the perturbative side has shown that observables like jet quenching [2] and elliptic flow [?] can also be described using a perturbative formalism. Since in phenomenological applications predictions are complicated by the modeling required to describe, for example, initial-state effects, the space-time evolution of the plasma, and hadronization of the plasma, there are significant theoretical uncertainties remaining. Therefore, one is hard put to conclude whether the plasma is strongly or weakly coupled based solely on RHIC data. To have a cleaner testing ground one can compare theoretical calculations with results from lattice quantum chromodynamics (QCD).
Looking forward to the upcoming heavy-ion experiments scheduled to take place at the Large Hadron Collider (LHC) at the European Organization for Nuclear Research (CERN)

JHEP00(2010)000
it is important to know if, at the higher temperatures generated, one expects a stronglycoupled (liquid) or weakly-coupled (plasma) description to be more appropriate. At RHIC, initial temperatures on the order of one to two times the QCD critical temperature, T c ∼ 190 MeV, were obtained. At LHC, initial temperatures on the order of 4−5 T c are expected.
The key question is, will the generated matter behave more like a plasma of quasiparticles at these higher temperatures.
The calculation of thermodynamic functions using weakly-coupled quantum field theory has a long history [4,5,6,7,8,9,10,11,12,13,14]. The QCD free energy is known up to order g 6 log(g); however, the resulting weak-coupling approximations do not converge at phenomenologically relevant couplings. For example, simply comparing the magnitude of low-order contributions to the QCD free energy with three quark flavors (N f = 3), one finds that the g 3 contribution is smaller than the g 2 contribution only for g ∼ 0.9 (α s ∼ 0.07) which corresponds to a temperature of T ∼ 10 5 GeV ∼ 5 × 10 5 T c .
In Fig 1, we show the weak-coupling expansion for the pressure of pure-glue QCD normalized to that of an ideal gas through order α 5/2 s . The various approximations oscillate wildy and show no signs of convergence in the temperature range shown. The bands are obtained by varying the renormalization scale µ by a factor of two around the value µ = 2πT and we use three-loop running of α s This oscillating behavior is generic for hot field theories and not specific to QCD. The poor convergence of finite-temperature perturbative expansions of thermodynamic functions stems from the fact that at high temperature the classical solution is not described by massless gluons. Instead one must include plasma effects such as the screening of electric fields and Landau damping of excitations via a self-consistent hard-thermal-loop (HTL) JHEP00(2010)000 resummation [15]. There are several ways of systematically reorganizing the perturbative expansion [16]. Here we will present the details of a new NNLO calculation which uses the hard-thermal-loop perturbation theory method [17,18,19,22] and compare with previously obtained LO and NLO results.
The basic idea of the technique is to add and subtract an effective mass term from the bare Lagrangian and to associate the added piece with the free part of the Lagrangian and the subtracted piece with the interactions [23,24]. However, in gauge theories, one cannot simply add and subtract a local mass term since this would violate gauge invariance [25,26,27]. Instead, one adds and subtracts an HTL improvement term which modifies the propagators and vertices self-consistently so that the reorganization is manifestly gauge invariant [25].
In this paper we discuss the calculation of thermodynamic functions of a gas of gluons at phenomenologically relevant temperatures using hard-thermal-loop perturbation theory. We present results at leading order (LO), next-to-leading order (NLO), and nextto-next-to-leading order (NNLO) and compare with available lattice data [28,29] for the thermodynamic functions of SU(3) Yang-Mills theory. The calculation is based on a reorganization of the theory around hard-thermal-loop (HTL) quasiparticles. Our results indicate that the lattice data at temperatures T ∼ 2 − 3 T c are consistent with the quasiparticle picture. This is a non-trivial result since, in this temperature regime, the QCD coupling constant is neither infinitesimally weak nor infinitely strong with g ∼ 2, or equivalently α s = g 2 /(4π) ∼ 0.3. Therefore, we have a crucial test of the quasiparticle picture in the intermediate coupling regime.

HTL perturbation theory
The Lagrangian density for pure-glue QCD in Minkowski space is 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π. 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

JHEP00(2010)000
The HTL improvement term is where the covariant derivative is D µ = ∂ µ − igA µ and y µ = (1,ŷ) is a light-like four-vector, and . . . y represents the average over the directions ofŷ. The term (1.4) has the form of the effective Lagrangian that would be induced by a rotationally invariant ensemble of charged sources with infinitely high momentum. The parameter m D can be identified with the Debye screening mass. HTLpt is defined by treating δ as a formal expansion parameter. The HTL perturbation expansion generates ultraviolet divergences. In perturbative QCD, renormalizability constrains the ultraviolet divergences to have a form that can be cancelled by the counterterm Lagrangian ∆L QCD . We will demonstrate that renormalized perturbation theory can be implemented by including a counterterm Lagrangian ∆L HTL among the interaction terms in (1.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 [17,18] that it was possible to renormalize the NLO order HTLpt prediction for the free energy of QCD using only a vacuum counterterm, a Debye mass counterterm, and a fermion mass counterterm. In this paper we will show that this is also possible at NNLO. In particular, the only new counterterm we need to introduce is for the coupling constant α s , which coincides with its perturbative value giving rise to the standard one-loop running.
We find that the counterterms necessary to renormalize HTLpt at NNLO are We note that the counterterm in Eq. (1.5) coincides with the perturbative one-loop running. 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 term in (1.4) are included to all orders but then systematically subtracted out at higher orders in perturbation theory by the δm 2 D terms in (1.4). If we set δ = 1, the Lagrangian (1.3) reduces to the QCD Lagrangian (1.1).
If the expansion in δ could be calculated to all orders the final result would not depend on m D when we set δ = 1. However, any truncation of the expansion in δ produces results that depend on m D . Some prescription is required to determine m D as a function of T and α s . We will discuss several prescriptions in Sec. VI.

Diagrams for the thermodynamic potential
In this section, we list the expressions for the diagrams that contribute to the thermodynamic potential through order δ 2 in HTL perturbation theory. The diagrams are shown JHEP00(2010)000 in Figs. 2, and 3. A key to the diagrams is given in Fig. 4. The expressions here will be given in Euclidean space; however, in Appendix A we present the HTLpt Feynman rules in Minkowski space.   The thermodynamic potential at leading order in HTL perturbation theory for QCD JHEP00(2010)000 is Here, F 1a+1b is the contribution from the gluon and ghost diagrams shown on the first line of Fig. 2 The transverse and longitudinal HTL propagators ∆ T (P ) and ∆ L (P ) are given in (A.49) and (A.50). The leading-order vacuum counterterm ∆ 0 E 0 was determined in Ref. [17]: 3) The thermodynamic potential at NLO in HTL perturbation theory can be written as where ∆ 1 E 0 and ∆ 1 m 2 D are the terms of order δ in the vacuum energy density and mass counterterms: The contributions from the two-loop diagrams with the three-gluon and four-gluon vertices are where R = Q − P . For the corresponding diagrams, see the second line of Fig. 2. The contribution from the ghost diagram is The contribution from the HTL gluon counterterm diagram with a single gluon selfenergy insertion is (2.10)

JHEP00(2010)000
The thermodynamic potential at NNLO in HTL perturbation theory can be written as where ∆ 1 α s , ∆ 2 m 2 D , and ∆ 2 E 0 are the terms of order δ 2 in the coupling constant, mass, and vacuum energy density counterterms: (2.14) The contributions from the three-loop diagrams are given by where S = −(P + Q + R) andΠ µν (P ) is the one-loop gluon self-energy defined by the yellow box in Fig. 4.

JHEP00(2010)000
The contributions from the two-loop diagrams with a single self-energy insertion are The two-loop diagrams with a subtracted vertex is The contribution from the HTL gluon counterterm diagram with two gluon self-energy insertions is (2.28)

Expansion in m D /T
In the papers [17,18,19,20,21], the free energy was reduced to scalar sum-integrals. It was clear that evaluating these scalar sum-integrals exactly was intractable and the sumintegrals were calculated approximately by expanding them in powers of m D /T . We will follow 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 is taken to be of order g. The free energy can be divided into contributions from hard and soft momenta. In the one-loop diagrams, the contributions are either hard (h) or soft (s), while at the twoloop 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.

Hard contribution
For hard momenta, the self-energies are suppressed by m D /T relative to the propagators, so we can expand in powers of Π T (P ) and Π L (P ).
For the one-loop graphs (1a) and (1b), we need to expand to second order in m 2 D :

Soft contribution
The soft contribution in the diagrams (1a + 1b) arises from the P 0 = 0 term in the sumintegral. At soft momentum P = (0, p), the HTL self-energy functions reduce to Π T (P ) = 0 and Π L (P ) = m 2 D . The transverse term vanishes in dimensional regularization because there is no momentum scale in the integral over p. Thus the soft contributions come from the longitudinal term only and read We have kept the order-ǫ in the m 2 D and m 3 D terms, respectively in Eqs. (3.1) and (3.2) since they contribute in the counterterms at next-to-leading order.

Hard contribution
The one-loop graph with a gluon self-energy insertion (2d) has an explicit factor of m 2 D and so we need only to expand the sum-integal to first order in m 2 D :

Soft contribution
The soft contribution from (2d) arises from the P 0 = 0 term in the sum-integral. Only the longitudinal part Π L (P ) of the self-energy contributes and reads (3.4) JHEP00(2010)000

(hh) contributions
For hard momenta, the self-energies are suppressed by m D /T relative to the propagators, so we can expand in powers of Π T and Π L . The two-loop contribution was calculated in Ref. [18] and reads (3.5) Using the expressions for the sum-integrals listed in Appendix B, we obtain

(hs) contributions
In the (hs) region, the momentum P is soft. The momenta Q and R are always hard. The function that multiplies the soft propagator ∆ T (0, p), ∆ L (0, p) or ∆ X (0, p) can be expanded in powers of the soft momentum p. In the case of ∆ T (0, p), the resulting integrals over p have no scale and they vanish in dimensional regularization. The integration measure p scales like m 3 D , the soft 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 . The two-loop contribution was calculated in Ref. [18] and reads In order to facilitate the calculations, it proves useful to isolate the terms that are specific to HTL perturbation theory. After integrating by parts and using the results from Zhai and Kastening [11], we can write Using the expressions for the integrals and sum-integrals in Appendices B and C, we obtain JHEP00(2010)000

(ss) contribution
The (ss) contribution to the free energy is given by a two-loop calculation in electrostatic QCD (EQCD) in three dimensions. This calculation was carried out in Ref. [9] by Braaten and Nieto. Alternatively, one can isolate the (ss) contributions from the two-loop diagrams which were calculated by Arnold and Zhai in Ref. [6]. In Ref. [18], this contribution was calculated and agrees with earlier results. One finds We have kept the order ǫ in terms Eqs. 3.3 Next-to-next-to-leading order

Hard contribution
The one-loop graph with two gluon self-energy insertions (3m) is proportional to m 4 D and so must be expanded to zeroth order in m 2 (3.11)

Soft contribution
The soft contribution from (3m) arises from the P 0 = 0 term in the sum-integral. Only the longitudinal part Π L (P ) of the self-energy contributes and reads . (3.12)

(hs) contributions
We also need the (hs) contribution from the diagrams (3h) − (3l). Again we calculate their contributions by expanding the two-loop diagrams (2a) − (2c) to first order in m 2 D . This yields This yields

(ss) contribution
The (ss) contribution from the two-loop diagrams with a single self-energy insertion or vertex correction can be obtained by expanding the two-loop result in powers of m 2 D . This yields We have verified this by explicitly calculating the relevant diagrams.

(hhh) contribution
If all the three loop momenta are hard, we can obtain the m D /T expansion simply by expanding in powers of m 2 D . To obtain the expansion through order g 5 , we can use bare propagators and vertices. The contributions from the three-loop diagrams were first calculated by Arnold and Zhai in Ref. [6], and later by Braaten and Nieto [9]. One finds

(hhs) contributions
All the diagrams except (3f) are infrared finite in the limit m D → 0. This implies that the g 5 contribution is given by using a dressed longitudinal propagator and bare vertices. The ring diagram (3f) is infrared divergent in that limit. The contribution through g 5 is obtained by expanding in powers of self-energies and vertices and one obtains (3.18)

(hss) contribution
For all 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 contribute as follows (3.19)

(sss) contribution
The (sss) contribution to the free energy is given by a three-loop calculation of the free energy of Electrostatic QCD in three dimensions. This calculation was performed in Ref. [9]. Alternatively, one can isolate the (sss) contributions from the diagrams listed in Ref. [6]. The result is (3.20) The expression for the integrals are given in Appendix C. Adding Eqs. (C.11)-(C.23), the final result is Note that all the poles in ǫ cancel.

The thermodynamic potential
In this section we present the final renormalized thermodynamic potential explicitly through order δ 2 , aka NNLO. The final NNLO expression is completely analytic; however, there are some numerically determined constants which remain in the final expressions at NLO.

Leading order
The leading order thermodynamic potential is given by the contribution from the diagrams (1a) and (1b).
JHEP00 (2010)000 where F ideal is the free energy of a gas of N 2 c − 1 massless spin-one bosons andm D andμ are dimensionless variables: The complete expression for the leading order thermodynamic potential is given by adding the leading vacuum energy counterterm (2.3) to Eq. (4.1): (4.5)

Next-to-leading order
The renormalization contributions at first order in δ are Using the results listed in Eqs. (2.5) and (2.6) the complete contribution from the counterterm at first order in δ is The NLO thermodynamic potential reads (4.8)

Next-to-next-to-leading order
The renormalization contributions at second order in δ are

JHEP00(2010)000
Using the results listed above, we obtain (4.10) Adding the NNLO counterterms (4.10) to the contributions from the various NNLO diagrams we obtain the renormalized NNLO thermodynamic potential. We note that at NNLO all numerically determined coefficients of order ǫ 0 drop out and we are left with a final result which is completely analytic. The resulting NNLO thermodynamic potential is Note that if we use the weak-coupling value for the Debye mass m 2 D = 4πN c α s T 2 /3, the NNLO HTLpt result (4.11) is guaranteed to reduce to the weak-coupling result through order g 5 and we have checked that this is the case.

Mass prescriptions
The mass parameter m D in HTLpt is completely arbitrary. To complete a calculation, it is necessary to specify m D as function of g and T . In this section we will discuss several prescriptions for the mass parameter.

Variational Debye mass
The variational mass is given by the solution to the equation

JHEP00(2010)000
At leading order in HTLpt, the only solution is the trivial solution, i.e. m D = 0. In that case, it is natural to chose the weak-coupling result for the Debye mass. This was done in Ref. [17].
At NLO, the resulting gap equation has a nontrivial solution, which is real for all values of the coupling. At NNLO, the solution to the gap equation (5.2) is plagued by imaginary parts for all values of the coupling. The problem with complex solutions seems to be generic since it has also been observed in screened perturbation theory [24] and QED [19,20,21]. In those cases, however, it was complex only for small values of the coupling.

Perturbative Debye mass
At leading order 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.
The next-to-leading order correction to the Debye mass is determined by resummation of one-loop diagrams with dressed vertices. Furthermore, since it suffices to take into account static modes in the loops, the HTL-corrections to the vertices also vanish. The result, however, turns out to be logarithmically infrared divergent, which reflects the sensitivity to the nonperturbative magnetic mass scale. The result was first obtained by Rebhan [43] and reads where m mag is the nonperturbative magnetic mass. We will not use this mass prescription since it involves the magnetic mass which would require input from e.g. lattice simulations.

BN mass parameter m 2 E
In the previous subsection, we saw that the Debye mass is sensitive to the nonperturbative magnetic mass which is of order g 2 T . In QED, the situation is much better. The Debye mass can be calculated order by order in e using resummed perturbation theory. The Debye mass then receives contributions from the scale T and eT . Effective field theory methods and dimensional reduction can be conveniently used to calculate separately the contributions to m D from the two scales in the problem. The contributions to m D and other physical quantities from the scale T can be calculated using bare propagators and vertices. The contributions from the soft scale can be calculated using an effective threedimensional field theory called electrostatic QED. The parameters of this effective theory are obtained by a matching procedure and encode the physics from the scale T . The effective field theory contains a massive field A 0 that up to normalization can be identified with the zeroth component of the gauge field in QED. The mass parameter m E of A 0 gives the contribution to the Debye mass from the hard scale T and can be written as a power series in e 2 .

JHEP00(2010)000
For nonabelian gauge theories, the corresponding effective three-dimensional theory is called electrostatic QCD. The mass parameter m E for the field A a 0 (which lives in the adjoint representation) can also be calculated as a power series in g 2 . It has been determined to order g 4 by Braaten and Nieto [9]. For pure-glue QCD, it reads We will use the mass parameter m E as another prescription for the Debye mass and denote it by the Braaten Nieto (BN) mass prescription.

Pressure
In this subsection, we present our results for the pressure using the variational mass prescription and the BN mass prescription. In Fig. 5, we compare the LO, NLO, and NNLO predictions for the real part of the pressure normalized to that of an ideal gas. Shaded bands show the result of varying the renormalization scale µ by a factor of two around µ = 2πT . We use three-loop running of α s .

Variational mass
In Fig. 6, we show the NNLO result for the imaginary part of the pressure normalized by the ideal gas pressure. We use three-loop running of α s . The imaginary part decreases with increasing temperature and is rather small beyond 3-4 T c . JHEP00(2010)000

BN mass
In Fig. 7, we show the HTLpt predictions for the pressure normalized to that of an ideal gas using the BN mass prescription. The bands are obtained by varying the renormalization JHEP00(2010)000 scale by a factor of two around µ = 2πT . We use one-loop running of α s . In Fig. 8, we again plot the normalized pressure, but now with three-loop running of α s . The agreement between the lattice data from Boyd et al. [28] is very good down to temperatures of around 3 T c . Comparing Figs. 7-8 we see that using the three-loop running, the band becomes wider. However, the difference is significant only for low T , where the HTLpt results disagrees with the lattice anyway. For T > 3T c , the prescription for the running makes very little difference. Until recently, lattice data for thermodynamic variables only existed for temperatures up to approximately 5 T c . In the paper by Enrodi et al [29], the authors calculate the pressure on the lattice for pure-glue QCD at very large temperatures. In Fig. 9, we show the results of Enrodi et al as well as Boyd et al, together with the HTLpt NLO and NNLO predictions for the pressure. The two points from Ref. [29] have large error bars, but the data points are consistent with the HTL predictions.
We note that our prediction for the normalized free energy using either the leadingorder, BN, or variational mass prescriptions is independent of N c if one holds α s N c fixed ('t Hooft limit). However, this need not be the case for an arbitrary mass prescription. The N c -independence of the normalized free energy of the free energy is in agreement with recent lattice measurements that show a very small dependence of the free energy on N c [30,31].
Due to the imaginary parts we do not show results coming from the variational prescription in the remainder of the results section. We will return to a discussion of the dependence on mass prescriptions in the conclusions. JHEP00(2010)000 Figure 9: Comparison of NLO, and NNLO predictions for the scaled pressure with SU(3) pure-glue lattice data from Boyd et al. [28] and Endrodi et al. [29]. Shaded bands show the result of varying the renormalization scale µ by a factor of two around µ = 2πT .

Energy density
The energy density E is defined by In Fig. 10, we show the LO, NLO, and NNLO predictions for energy density normalized to that of an ideal gas. We use three-loop running and the BN mass. The bands show the result of varying the renormalization scale µ by a factor of two around µ = 2πT . Our NNLO predictions are in excellent agreement with the lattice data down to T ≃ 2T c .

Entropy
The entropy density is defined by where all other parameters that F depends on, are kept fixed. In Fig. 11, we show the entropy density normalized to that of an ideal gas. We use three-loop running and the BN mass. The points are lattice data from Boyd et al. [28]. Our NNLO predictions are in excellent agreement with the lattice data down to T ≃ 2T c .

Trace anomaly
In pure-glue QCD or in QCD with massless quarks, there is no mass scale in the Lagrangian and the theory is scale invariant. At the classical level, this implies that the trace of the energy-momentum tensor vanishes. At the quantum level, scale invariance is broken by JHEP00(2010)000  be written as In Fig. 12, we show the HTLpt predictions for the trace anomaly divided by E ideal using the BN mass prescription and three-loop running of α s . The points are lattice data from Boyd et al. [28] For temperatures below approximately 2T c , there is a large discrepancy between the HTLpt predictions and lattice data. At LO and NLO, the curves are even bending downwards. At temperatures close to the phase transition it has been suggested that the discrepancy between HTLpt resummed predictions for thermodynamics functions and, in particular, the trace anomaly is due to influence of a power corrections [32,33,34,35,36] which are related to confinement. Phenomenological fits of lattice data which include such power corrections show that the agreement with lattice data is improved [37,38]. Alternatively, others have constructed AdS/CFT inspired models which break conformal invariance "by hand" [39,40,41,42]. These models are also able to fit the thermodynamical functions of QCD at temperatures close to the phase transition. Figure 12: Comparison of LO, NLO, and NNLO predictions for the scaled trace anomaly with SU(3) pure-glue lattice data from Boyd et al. [28]. Shaded bands show the result of varying the renormalization scale µ by a factor of two around µ = 2πT .
In Fig. 13, we show the HTLpt predictions for the trace anomaly scaled by T 2 /T 2 c using the BN mass prescription and three-loop running of α s . The points are lattice data from Boyd et al. [28]. The most remarkable feature is that lattice data are essentially constant over a very large temperature range. Clearly, HTLpt does not reproduce the scaled lattice data precisely; however, the agreement is dramatically improved when going from NLO to NNLO. JHEP00(2010)000 Figure 13: Comparison of LO, NLO, and NNLO predictions for the scaled trace anomaly multiplied by T 2 /T 2 c with SU(3) pure-glue lattice data from Boyd et al. [28]. Shaded bands show the result of varying the renormalization scale µ by a factor of two around µ = 2πT . Orange band shown with lattice data indicates scaled error assuming a ±2.5% error in the lattice data for the trace anomaly.

Conclusions
In this paper, we have presented results for the LO, NLO, and NNLO thermodynamic functions for SU(N c ) Yang-Mills theory using HTLpt. We compared our predictions with lattice data for N c = 3 and found that with a perturbative mass prescription HTLpt is consistent with available lattice data down to approximately T ∼ 3 T c in the case of the pressure and T ∼ 2 T c in the case of the energy density and entropy. These results are in line with expectations since below T ∼ 2 − 3 T c a simple "electric" quasiparticle approximation breaks down due to nonperturbative chromomagnetic effects.
The mass parameter m D in HTLpt is arbitrary and we employed two different prescriptions for fixing it. Unfortunately, the variational gap equation has four complex conjugate solutions, two with positive real parts. This has also been observed in scalar theory and QED. Whether this is a problem of HTLpt as such or is related to our m D /T expansion is unknown. Since it is not currently possible to evaluate the NNLO HTLpt diagrams in gauge theories exactly, it is impossible to settle the issue at this stage. On the other hand, the BN mass prescription is well defined to all orders in perturbation theory and does a reasonable job reproducing available lattice data for temperatures above T ≃ 3T c .
That being said without lattice data to compare with one would be hard pressed to favor one prescription over the other. While it is true that variational solutions are complex, one could be tempted to ignore the imaginary contributions and take Re[P] as the approximation. As we have shown (see Fig. 5) in this case the convergence of HTLpt is impressive; however, the result lies above the lattice data at low temperatures. Taking JHEP00(2010)000 the BN mass prescription on the other hand results in a real-valued pressure which seems to be in better agreement with the lattice data at the expense of poorer convergence of the successive approximations. The dependence on the mass prescription gives another measure of our theoretical uncertainty. If one were able to carry the HTLpt δ expansion to higher and higher orders the dependence on the mass prescription would go away; however, going to higher orders presents a rather daunting task since one encounters a purely nonperturbative contribution at O(α 3 s ) for which lattice input is required. Finally, we emphasize that HTLpt is gauge invariant by construction and that it can be used to calculate time-dependent quantities as well since it is formulated in Minkowski space. The fact that our BN mass results are consistent with lattice data down to approximately 3T c , suggests that HTLpt can provide a coherent framework that can be used to systematically calculate real-time quantities such as heavy quark diffusion and viscosity at temperatures relevant for LHC.

Acknowledgments
We thank H. Stöcker for his encouragement and support for this endeavor. N. Su

A. HTL Feynman rules
In this appendix, we present Feynman rules for HTL perturbation theory in QCD. We give explicit expressions for the propagators and for the 3-and 4-gluon vertices. The Feynman rules are given in Minkowski space to facilitate future application to real-time processes. A Minkowski momentum is denoted p = (p 0 , p), and the inner product is p · q = p 0 q 0 − p · q. The vector that specifies the thermal rest frame is n = (1, 0).

A.1 Gluon self-energy
The HTL gluon self-energy tensor for a gluon of momentum p is The tensor T µν (p, q), which is defined only for momenta that satisfy p + q = 0, is The angular brackets indicate averaging over the spatial directions of the light-like vector y = (1,ŷ). The tensor T µν is symmetric in µ and ν and satisfies the "Ward identity"

A.4 Four-gluon vertex
The four-gluon vertex for gluons with outgoing momenta p, q, r, and s, Lorentz indices µ, ν, λ, and σ, and color indices a, b, c, and d is where the cyclic permutations are of (q, ν, b), (r, λ, c), and (s, σ, d). The matrices T a are the fundamental representation of the SU (3) algebra with the standard normalization tr(T a T b ) = 1 2 δ ab . The tensor T µνλσ in the HTL correction term is defined only for p + q + r + s = 0: T µνλσ (p, q, r, s) = y µ y ν y λ y σ p·n p·y q·y (q + r)·y + (p + q)·n q·y r·y (r + s)·y + (p + q + r)·n r·y s·y (s + p)·y .
This tensor is totally symmetric in its four indices and traceless in any pair of indices: g µν T µνλσ = 0. It is even under cyclic or anti-cyclic permutations of the momenta p, q, r, and s. It satisfies the "Ward identity" q µ T µνλσ (p, q, r, s) = T νλσ (p + q, r, s) −T νλσ (p, r + q, s) (A. 38) and the "Bianchi identity" T µνλσ (p, q, r, s) + T µνλσ (p, r, s, q) + T µνλσ (p, s, q, r) = 0 . (A.39) When its color indices are traced in pairs, the four-gluon vertex becomes particularly simple: where the color-traced four-gluon vertex tensor is Γ µν,λσ (p, q, r, s) = 2g µν g λσ − g µλ g νσ − g µσ g νλ − m 2 D T µνλσ (p, s, q, r) . Every closed ghost loop requires a multiplicative factor of −1.

A.6 HTL counterterm
The Feynman rule for the insertion of an HTL counterterm into a gluon propagator is where Π µν (p) is the HTL gluon self-energy tensor given in (A.8).

A.7 Imaginary-time formalism
In the imaginary-time formalism, Minkoswski energies have discrete imaginary values p 0 = i(2πnT ) and integrals over Minkowski space are replaced by sum-integrals over Euclidean vectors (2πnT, p). We will use the notation P = (P 0 , p) for Euclidean momenta. The magnitude of the spatial momentum will be denoted p = |p|, and should not be confused with a Minkowski vector. 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). The Feynman rules for Minkowski space given above can be easily adapted to Euclidean space. The Euclidean tensor in a given Feynman rule is obtained from the corresponding Minkowski tensor with raised indices by replacing each Minkowski energy p 0 by iP 0 , where P 0 is the corresponding Euclidean energy, and multipying by −i for every 0 index. This prescription transforms p = (p 0 , p) into P = (P 0 , p), g µν into −δ µν , and p·q into −P·Q. The effect on the HTL tensors defined in (A.2), (A.33), and (A.37) is equivalent to substituting p·n → −P ·N where N = (−i, 0), p·y → −P ·Y where Y = (−i,ŷ), and y µ → Y µ . For example, the Euclidean tensor corresponding to (A.2) is The average is taken over the directions of the unit vectorŷ. Alternatively, one can calculate a diagram by using the Feynman rules for Minkowski momenta, reducing the expressions for diagrams to scalars, and then make the appropriate JHEP00(2010)000 substitutions, such as p 2 → −P 2 , p · q → −P · Q, and n · p → in · P . For example, the propagator functions (A.21) and (A.22) become . (A.50) The expressions for the HTL self-energy functions Π T (P ) and Π L (P ) are given by (A.13) and (A.14) with n 2 p replaced by n 2 P = p 2 /P 2 and T 00 (p, −p) replaced by Note that this function differs by a sign from the 00 component of the Euclidean tensor corresponding to (A.2): T 00 (P, −P ) = −T 00 (p, −p) A more convenient form for calculating sum-integrals that involve the function T P is where the angular brackets represent an average over c defined by and w(ǫ) is given in (A.16).

B. Sum-integrals
In the imaginary-time formalism for thermal field theory, the 4-momentum P = (P 0 , p) is Euclidean with P 2 = P 2 0 + p 2 . The Euclidean energy P 0 has discrete values: P 0 = 2nπT for bosons, where n is an integer. Loop diagrams involve sums over P 0 and integrals over p. With dimensional regularization, the integral is generalized to d = 3−2ǫ spatial dimensions. We define the dimensionally regularized sum-integral by where 3 − 2ǫ is the dimension of space and µ 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. Below we list the sum-integrals required to complete the three loop calculation. We refer to Ref. [18] for details concerning the sum-integral evaluations.

B.1 One-loop sum-integrals
The simple one-loop sum-integrals required in our calculations can be derived from the formulas (B. 2) The specific bosonic one-loop sum-integrals needed are P log P 2 = − π 2 45 T 4 , (B.3) The number γ 1 is the first Stieltjes gamma constant defined by the equation We also need some more difficult one-loop sum-integrals that involve the HTL function defined in (A.51). The specific bosonic sum-integrals needed are where R = −(P + Q) and r = |p + q|. We also need some more difficult two-loop sumintegrals that involve the functions T P defined in (A.51). The specific bosonic sum-integrals needed are P Q 1 P 2 Q 2 r 2 T R = The three-loop sum-integrals were first calculated by Arnold and Zhai and calculational details can be found in Ref. [6].

C. Three-dimensional integrals
Dimensional regularization can be used to regularize both the ultraviolet divergences and infrared divergences in 3-dimensional integrals over momenta. The spatial dimension is generalized to d = 3 − 2ǫ dimensions. Integrals are evaluated at a value of d for which they converge and then analytically continued to d = 3. We use the integration measure p ≡ e γ E µ 2 4π ǫ d 3−2ǫ p (2π) 3−2ǫ .

C.1 One-loop integrals
The one-loop integral is given by Specifically, we need (C.5)

C.2 Two-loop integrals
We also need a few two-loop integrals on the form Specifically, we need J 1 , J 2 , and K 1 which were calculated in Refs. [9,11]: (C.10)

C.3 Three-loop integrals
We also need a number of three-loop integrals. The specific integrals we need are listed below and were calculated in Refs. [9,11]. They are special cases of more general integrals defined in Ref. [44].