Two-loop HTL-resummed thermodynamics for N=4 supersymmetric Yang-Mills theory

We compute the two-loop hard-thermal-loop (HTL) resummed thermodynamic potential for N=4 supersymmetric Yang-Mills (SYM). Our final result is manifestly gauge-invariant and was renormalized using only simple vacuum energy, gluon mass, scalar mass, and quark mass counter terms. The HTL mass parameters m_D, M_D, and m_q are then determined self-consistently using a variational prescription which results in a set of coupled gap equations. Based on this, we obtain the two-loop HTL-resummed thermodynamic functions of N=4 SYM. We compare our final result with known results obtained in the weak- and strong-coupling limits. We also compare to previously obtained approximately self-consistent HTL resummations and Pad\'e approximants. We find that the two-loop HTL resummed results for the scaled entropy density is a quantitatively reliable approximation to the scaled entropy density for 0<= lambda<~ 2 and is in agreement with previous approximately self-consistent HTL resummation results for lambda<~ 6.

1 Introduction N = 4 supersymmetric Yang-Mills theory (SYM) is the most famous example of a conformal field theory (CFT) in four dimensions, and is often taken as a model for hot QCD in the large number of colors N c and strong 't Hooft coupling λ = g 2 N c limits. The strong coupling behavior of the free energy has been computed using the anti-de Sitter space/CFT (AdS/CFT) correspondence [1], with the result being where F ideal = −d A π 2 T 4 /6 is the ideal or Stefan-Boltzmann limit of the free energy and S ideal = 2d A π 2 T 3 /3, with d A = N 2 c − 1 being the dimension of the adjoint representation. In the weak-couping limit the N = 4 SYM free energy has been calculated through order λ 3/2 giving [2][3][4] (1.2) Note that, since the beta function of the N = 4 SYM theory is zero, the coupling constant does not run and is independent of the temperature. As a result, we can vary the coupling between the two limits at each temperature. One expects these two series to describe their respective asymptotic limits correctly, however, the radius of convergence of each of these series is unknown and, therefore, it Weak and strong coupling results for the entropy density in N = 4 SYM theory compared to a R [4,4] Padè approximation constructed from both limits.
is unclear to what degree each of these can trusted away from their respective limits. In Fig. 1 we plot the scaled entropy density resulting from Eqs. (1.1) and (1.2) as a function of λ along with a R [4,4] Padé approximant constructed from these results [3,5]. 1 From this Figure, we can see that the two successive weak coupling approximations are only close to one another below λ ∼ 0.1 and rapidly diverge beyond λ ∼ 1. In the strong coupling limit, only the first two terms in the series are known. As can be seen from Fig. 1 the strong coupling result diverges quickly below λ ∼ 10. The question then becomes, how can we systematically extend these two results into the intermediate coupling region λ ∼ 1 − 10. In this paper, we present progress towards this goal in the weak-coupling limit using hard-thermal-loop (HTL) perturbation theory. Our study is complementary to the earlier work of Blaizot, Iancu, and Rebhan in which they applied HTL resummation using an approximately self-consistent scheme [5]. The key difference from this earlier work is that our result is manifestly gauge invariant and based on the systematically improvable HTL perturbation theory (HTLpt) framework [6][7][8].
The goal of our work is to improve the convergence of the successive weak-coupling approximations. One promising approach is to use a variational framework in which the free energy F is expressed as the variational minimum of the thermodynamic potential Ω(T, λ; m 2 ) that depends on one or more variational parameters that we denote collectively by a F(T, λ) = Ω(T, λ; a) ∂Ω/∂a=0 . (1.3) For example, the Φ-derivable approximation is a widely used variational method in which 1 Although a Padé approximant might provide a convenient interpolation between the weak-and strongcoupling limits their construction is in no sense systematic. In particular, the resulting expressions are incomplete since we know that, at least at weak coupling, the series will contain logarithms of the coupling constant beyond O(λ 3/2 ). the propagator is used as an infinite set of variational parameters [9,10]. The Φ-derivable thermodynamic potential Ω is given by the 2-particle-irreducible (2PI) effective action, which is the sum of all diagrams that are 2-particle-irreducible with respect to the complete propagator [11]. This method is difficult to apply to relativistic field theories except for the case where the self-energy is momentum-independent. Despite this, there still have some progress in applications to quantum chromodynamics (QCD) and electrodynamics (QED) [12][13][14][15][16][17]. Historically, the Φ-derivable approximation was first applied to QCD by Freedman and McLerran [18], who demonstrated that the thermodynamic potential Ω is gauge dependent beyond a given order in the coupling constant. The gauge parameter dependence appears at the same order in α s as the series truncation when evaluated off the stationary point and at twice the order in α s when evaluated at the stationary point [17,19,20]. Despite this issue, this method had been used as the starting point for approximately self-consistent HTL resummation of the entropy [21,22] and the pressure [23]. The problems encountered when applying the Φ-derivable approximation to gauge theories motivated the use of alternative variational approximations. One such alternative, which in its simplest form involves a single variational parameter m, has been called optimized perturbation theory [24], variational perturbation theory [25,26], or the linear δ expansion [27,28]. This strategy has been successfully used for the thermodynamics of the massless φ 4 field theory up to four-loop order using "screened perturbation theory" [29][30][31][32], and spontaneously broken field theories at finite temperature [33][34][35][36]. One impediment to applying such ideas to gauge theories was that one cannot simply introduce a scalar mass for the gluon without breaking gauge invariance. The solution to this problem was introduced in Refs. [7] and [8] in which it was demonstrated that one could generalize the linear delta expansion by adding and subtracting the full gauge-invariant HTL effective Lagrangian [37,38].
The resulting scheme was called hard-thermal-loop perturbation theory (HTLpt) [7,8]. The HTLpt method has been used to improve the convergence of weak coupling calculations of the free energy in QED [39] and QCD up to three-loop order at finite temperature and chemical potential [40][41][42][43][44][45][46]. When confronted with finite temperature and chemical potential lattice QCD data the HTLpt resummation scheme works remarkably well down to temperatures on the order of 200-300 MeV where the QCD coupling constant is on the order of g s ∼ 2 [45][46][47]. The method successfully describes all thermodynamic variables including second-and fourth-order quark susceptibilities. Herein, we will take the first steps in applying this method to N = 4 SYM in the hope that a similar improvement in convergence can be achieved in this theory at intermediate couplings. We will calculate the one-and two-loop HTLpt-resummed thermodynamic potential Additionally, in N = 4 SYM using the same method as was used to obtain the one-and two-loop QCD results in Refs. [7,8,48,49]. Importantly, in these papers it was demonstrated that it was possible to renormalize the resummed thermodynamic potential at two-loop order using only known vacuum and mass counterterms. Herein we will demonstrate the same occurs in N = 4 SYM. In this theory the NLO contributions include scalar and scalar-gluon, scalar-quark interactions compared to the QCD calculation, however these are relatively straightforward to include. In N = 4 SYM instead of having only gluon and quark thermal masses, m D and m q , respectively, we will also have a thermal mass for the scalar particles, M D . Our results at one-and two-loop order are infinite series in λ which, when trucated at O(λ 3/2 ), reproduce the weak-coupling results obtained previously in Refs. [2][3][4]. In order to make the calculation tractable, we expand the HTLpt scalarized sum-integrals in a power series in the three mass parameters M D , m D , and m q such that it includes terms that would naively contribute throughx O(λ 5/2 ). Our final results indicate that, in N = 4 SYM, NLO HTLpt provides a good approximation for the scaled entropy for couplings in the range 0 ≤ λ 2.
We begin with a brief summary of HTLpt for N = 4 SYM in Sec. 2. In Sec. 3, we give the expressions for the one-and two-loop diagrams contributing to the SYM thermodynamic potential. In Sec. 4, we reduce these diagrams to scalar sum-integrals. As mentioned in the prior paragraph, since it would be intractable to calculate the resulting sum-integrals analytically, in Sec. 5 we expand these expressions by treating m D , M D and m q as O(λ 1/2 ) and expanding the integrals in powers of m D /T , M D /T and m q /T , keeping all terms up to O(λ 5/2 ). In Sec. 6, we combine the results obtained in Sec. 5 to obtain the complete expressions for the leading-(LO) and next-to-leading order (NLO) thermodynamic potentials. In Sec. 7, we present our numerical results for the HTLpt-resummed LO and NLO scaled thermodynamic functions in N = 4 SYM and compare to prior results in the literature. For details concerning the transformation to Euclidean space and the sum-integrals necessary we refer the reader to the appendices of Refs. [48,49].

Notation and conventions
We use lower-case letters for Minkowski space four-vectors, e.g. p, and upper-case letters for Euclidean space four-vectors, e.g. P . We use the mostly minus convention for the metric.

HTLpt for N = 4 SYM
In N = 4 SYM theory all fields belong to the adjoint representation of the SU (N c ) gauge group. For the fermionic fields, a massless two-component Weyl fermion ψ in four dimensions can be converted to four-component Majorana fermions [50][51][52][53][54] where α = 1, 2 and the Weyl spinors satisfyψα ≡ [ψ α ] † . The conjugate spinorψ is not independent, but is related to ψ via the Majorana condition ψ = Cψ, where C = αβ 0 0 αβ is the charge conjugation operator with 02 = − 11 ≡ −1. In the following, we will use the indices i, j = 1, 2, 3, 4 to enumerate the Majorana fermions and use ψ i to denote each bispinor. The definition of gauge field is the same as QCD, and A µ can be expanded as A µ = A a µ t a , with real coefficients A a µ , and Hermitian color generators t a in the fundamental representation that satisfy where a, b = 1, · · · , N 2 c − 1, the structure constants f abc are real and completely antisymmetric. The fermionic fields can similarly be expanded in the basis of color generators as ψ i = ψ a i t a . The coefficients ψ a i are four-component Grassmann-valued spinors. There are six independent real scalar fields which are represented by a multiplet where X p and Y q hermitian, with p,q = 1, 2, 3. X p and Y q denote scalars and pseudoscalar fields, respectively. We will use a capital Latin index A to denote components of vector Φ. Therefore Φ A , X p , and Y q can be expanded as Φ A = Φ a A t a , with A = 1, · · · , 6, and The Lagrangian density that generates the perturbative expansion for N = 4 SYM theory in Minkowski-space can be expressed as where the field strength tensor is is the covariant derivation in the adjoint representation. α p and β q are 4 × 4 matrices that satisfy and their explicit form can be given as where σ i with i ∈ {1, 2, 3} are the 2×2 Pauli matrices. And α and β satisfies α p ik α p kj = −3δ ij and β q ij β p ji = −4δ pq , with δ ii = 4 for four Majorana fermions and δ pp = 3 for three scalars. The ghost term L gh depends on the choice of the gauge-fixing term L gf and is the same as in QCD. Here we work in general covariant gauge, giving with ξ being the gauge parameter.
In general, perturbative expansion in powers of g in quantum field theory generates ultraviolet divergences. The renormalizability of perturbation theory guarantees that all divergences in physical quantities can be removed by the renormalization of masses and coupling constants. The coupling constant in N = 4 SYM theory is denoted as λ = g 2 N c . Unlike QCD, N = 4 SYM theory does not run.
Similar to the case of QCD presented in Refs. [48,49], HTLpt is also a reorganziation of the perturbation series for the SYM theory, and can be defined by introducing an expansion parameter δ. The HTL "shifted" Lagrangian density can be written as (2.8) The HTL improvement term is where y µ = (1,ŷ) is a light-like four vector defined in App. A, j ∈ {1 . . . 4} indexes the four Majorana fermions, A ∈ {1 . . . 6} indexes the scalar degrees of freedom, and · · · ŷ represents the average over the direction ofŷ. The parameters m D and M D are the electric screening masses for the gauge field and the adjoint scalar field, respectively. The parameter m q can be seen as the induced finite temperature quark mass. We note that, if we set δ = 1, the Lagrangian above (2.8) reduces to the vacuum N = 4 SYM Lagrangian (2.4). HTLpt is defined by treating δ as a formal expansion parameter, expanding around δ = 0 to a fixed order, and then setting δ = 1. In the limit that this expansion is taken to all orders, one reproduces the QCD result by construction, however, the loop expansion is now shifted to be around the high-temperature minimum of the effective action, resulting in a reorganization of the perturbation series which has better convergence than the naive expansion loop expansion around the T = 0 vacuum. In addition, this reorganization eliminates all infrared divergences associated with the electric sector of the theory. The HTLpt reorganization generates new ultraviolet (UV) divergences and, due to the renormalizability of perturbation theory, the ultraviolet divergences are constrained to have a form that can be canceled by the counterterm Lagrangian ∆L HTL . References [48,49] demonstrated that at two-loop order the thermodynamic potential can be renormalized using a simple counterterm Lagrangian ∆L HTL containing vacuum and mass counterterms. Although the general structure of the ultraviolet divergences is unknown, it has been demonstrated that one can renormalize the next-to-leading order HTLpt thermodynamic potential through three-loop order using only vacuum, gluon thermal mass, quark thermal mass, and coupling constant counterterms [43,46]. In this paper, we demonstrate that the same method can be used for N = 4 SYM and we compute the vacuum and screening mass counterterms necessary.
We find that the vacuum counterterm ∆ 0 E 0 , which is the leading order counterterm in the δ expansion of the vacuum energy E 0 , can be obtained by calculating the free energy to leading order in δ. In Sec. 6.1, we show that ∆ 1 E 0 can be obtained by expanding ∆E 0 to linear order in δ. As a result, the counterterm ∆E 0 has the form To calculate the NLO free energy we need to expand to order δ and we will need the counterterms ∆E 0 , ∆m 2 D , ∆m 2 q , and ∆M 2 D to order δ in order to cancel the UV divergences. We find that in order to remove the divergences to two-loop order, the mass counterterms should have the form In the N = 4 SYM theory, we will use the same method as in QCD to calculate physical observables in HTLpt, namely expanding the path-integral in powers of δ, truncating at some specified order, and then setting δ = 1. The results of the physical observables will depend on m D , M D , and m q for any truncation of the expansion in δ, and some prescription is required to determine m D , M D , and m q as a function of λ. In this work, we will follow the two-loop HTLpt QCD prescription and determine them by minimizing the free energy. If we use Ω N (T, λ, m D , M D , m q , δ) to represent the thermodynamic potential expanded to N -th order in δ, then our full variational prescription is We will call Eqs. (2.12) the gap equations. The free energy is obtained by evaluating the thermodynamic potential at the solution to the gap equations. Other thermodynamic functions can then be obtained by taking appropriate derivatives of free energy with respect to T .

Next-to-leading order thermodynamic potential
In the imaginary-time formalism, Minkoswski energies have discrete imaginary values p 0 = i(2πnT ), and the integrals over Minkowski space should be replaced by Euclidean sum integrals. There are two ways to do this which have been discussed in Refs. [48,49]. One is transforming the Feynman rules in Minkowski space given in App. A into the form in Euclidean space firstly, then calculating the free energy. The other way is using the Feynman rules in Minkowski space to get the forms of free energy, after reducing these forms, transforming it into the form in Euclidean space. Results from the two methods must be the same. The HTL perturbative thermodynamic potential at next-to-leading order in N = 4 SYM can be expressed as where Ω LO is the leading order thermodynamic potential, O(δ 0 ), which includes the oneloop graphs shown in Fig. 2 and the LO vacuum renormalization counterterm. We first discuss the contributions at this order.

LO thermodynamic potential
In SU(N c ) gauge theory with massless particles, Ω LO can be expressed as There are four independent Majorana fermions in the adjoint representation, d F = 4d A , and d S = 6d A for the six scalars.
There are D = d + 1 polarization state for gluons, where d is the number of spatial dimensions. After canceling the two unphysical states using the ghost contribution, we obtain the HTL one-loop free energy of each of the color states of the gluon The transverse and longitudinal HTL propagators ∆ T (P ) and ∆ L (P ) are the HTL gluon propagator (A.13) in Euclidean space , .
The result above is the same as in QCD. The only difference is the definition of m 2 D in the gluon propagator which now contains contributions from gluon, fermion, and scalar loops as detailed in Appendix A.1.
Since the Majorana fermion is its own antiparticle, the fermionic contribution is reduced by a factor of two when comparing QCD and N = 4 SYM. Our definition of m 2 q is presented in Appendix A.3. The one-loop fermionic free energy is where Σ(P ) is the HTL fermion self-energy (A.27) in Euclidean space. The scalar one-loop free energy is simply

5)
where ∆ −1 s (P ) is the inverse scalar propagator which is given in Eq. (4.2). Finally, we note that the leading order counterterm ∆ 0 E 0 cancels the divergent terms of the one-loop thermodynamic potential in N = 4 SYM theory.

NLO thermodynamic potential
In Eq. (3.1) Ω 2-loop corresponds to the two-loop contributions shown in Fig. 3. It can be expressed as where λ = g 2 N c is the 't Hooft coupling constant.
The gluon propagator, the three-gluon vertex, the four-gluon vertex, and the gluonghost vertex are the same as in QCD up to the expression for the Debye mass m D . 2 As a result, the purely gluonic and glue-ghost graphs given by F 3g , F 4g , and F gh are, respec-tively, where R + P + Q = 0. Contributions come from the two-loop diagrams with four scalar vertex, scalar-gluon vertex and scalar-gluon four vertex are respectively where R + P − Q = 0. The contributions F 3qg and F 4qg involve only quarks and gluons and, since the Majorana fermion is its own antiparticle, their symmetry factor is 1/4 instead of 1/2 in QCD. Additionally, there are four Majorana fermions in N = 4 SYM, so that F 3qg and F 4qg are 2 times the result obtained in QCD. As a result, we can substitute s F and d F in Ref. [49] to 2N c and 2d A , respectively, to obtain the N = 4 SYM result. After this adjustment, the only other change required is to use the N = 4 SYM definitions of m 2 D and m 2 q . Based on the results contained in Ref. [49] one obtains The momentum conservation is R + P − Q = 0. One can also obtain these expressions using the Feynman rules contained in App. A. The final new graph, the quark-scalar diagram F 3qs , can be split into two parts, one coming from the quark-scalar vertex, and the other coming from the quark-pseudoscalar vertex. Using the Feynman rules in App. A, one finds that their contributions are the same. As a consequence, F 3qs can be written as where R + P − Q = 0.
The contribution Ω HTL in Eq. (3.1) is the sum of the gluon, quark and scalar HTL counterterms shown in Fig. 3. These enter in order to subtract contributions at lower loop orders and guarantee that naive perturbative results are recovered order by order if the expressions are truncated in λ. They can be expressed as (3.12) There are two ways to get these three contributions, one is using the Feynman rules in Appendix A, the other one is substituting in the one-loop expressions for F g +F ghost , F s , and F q and expanding them to linear order in δ. In terms of the first method, the contribution from the HTL gluon counterterm diagram is It is the same as in QCD up to the definition of m 2 D . The contribution from the HTL scalar counterterm diagram is where P AA aa (P ) = M 2 D is referred in Appendix A.6. The contribution from the HTL quark counterterm diagram is Compared to QCD this is different by one half due to the Majorana nature of the SYM fermions. As usual, the quark mass should be adjusted to the SYM case.
Since HTL perturbation theory is renormalizable, the ultraviolet divergences of free energy at any order in δ can be cancelled by E, m 2 D , m 2 q , and M 2 D and the coupling constant λ. ∆Ω NLO in Eq. (3.1) is the renormalization contribution at first order in δ is used to cancel the next-to-leading order divergences. It can be expressed as where ∆ 1 E 0 , ∆ 1 m 2 D , ∆ 1 m 2 q , and ∆ 1 M 2 D are the terms of order δ in the vacuum energy (2.10) and mass counterterms (2.11). The first term ∆ 1 E 0 can be obtained simply by expanding ∆ 0 E 0 to first order in δ. The second term in (3.16) is slightly different from the QCD result in Refs. [48,49]. This is because this term must be used to cancel the divergences of two-loop self energy. As we can see in (6.5), there are two mixed term M D m 2 D and M D m 2 q which comes from the (hs) contribution of F 4gs and F 3qs respectively. There are many ways to construct the mass renormalization form which is corresponding to the second part of eq.(3.16), but the only way to use one set of three counterterms ∆ 1 m 2 D , ∆ 1 m 2 q , and ∆ 1 M 2 D is the form we have shown above.
In this work, we calculate the thermodynamic potential as an expansion in powers of m D /T , m q /T , and M D /T to order g 5 . We will show that, at order δ, all divergences in the two-loop thermodynamic potential plus HTL counterterms can be removed by these vacuum and mass counterterms. This means the method used in QCD can also be used in N = 4 SYM theory. This provides nontrivial evidence for the renormalizability of HTLpt at order δ in N = 4 SYM.

Reduction to scalar sum-integrals
Since we can make use of prior QCD results, we only need to calculate F s , F sct , F 4s , F 3gs , F 4gs , F 3qs , and the HTL counterterms contributing to Eq. (3.12). The first step to calculate the new SYM contributions in Figs. 2 and 3 is to reduce the sum of these diagrams to scalar sum-integrals. In Euclidean space, by substituting p 0 to iP 0 the scalar propagator can be written as so its inverse is The leading-order scalar contribution can be written as The HTL scalar counterterm can be written as We proceed to simplify the sum of formulas in Eq. (3.9) in general covariant gauge parameterized by ξ. Using Eqs. (4.1) and (A.24), we obtain F 4s+3gs+4gs = 15 2 P Q ∆ s (P )∆ s (Q) There are two terms which depend on ξ in (4.5), however, using P + R − Q = 0, one finds that they cancel each other, so that the sum of these contributions is gauge independent. Similarly, the results for F 3g+4g+gh and F 3qg+4qg are independent of gauge parameter ξ as shown in prior QCD calculations. Therefore, we have verified explicitly that the NLO HTLpt resummed thermodynamic potential in N = 4 SYM is gauge independent. Finally, we simplify Eq. (3.11) to where in Euclidean space with T P is the angular average T 00 (p, −p) in Euclidean space, can be expressed as where the angular brackets denote an average over c defined by where ω( ) is given in (A.15).

High temperature expansion
Having reduced F s , F 4s+3gs+4gs , F 3qs , and the HTL counterterm F sct to scalar sumintegrals, we will now evaluate these sum-integrals approximately by expanding them in powers of m D /T , m q /T , and M D /T . We will keep terms that contribute through O(g 5 ) if m D , m q and M D are taken to be of order g at leading-order. Additionally, each of these terms can be divided into contributions from hard and soft momentum, so we will proceed to calculate their hard and soft contributions, respectively. In some cases, the results presented here were obtained in previous one-and two-loop QCD HTLpt papers [8,48,49]. When converting the prior QCD graphs involving quarks, as mentioned previously, one has to take into account that the four SYM quarks are Majorana fermions. Here we list results for all contributions to the N = 4 SYM Feynman graphs and counterterms for completeness and ease of reference. In all cases, we use the integral and sum-integral formulas from Refs. [48,49] to obtain explicit expressions.

One-loop sum-integrals
The one-loop sum-integrals include the leading gluon, quark, and scalar contributions (3.3), (3.5), and (3.6) along with their corresponding counterterms (3.13), (3.15), and (3.14). In order to include all terms through O(g 5 ), we need to expand the one-loop contribution to order m 4 D , m 4 q , and M 4 D .

Hard contributions
The hard contribution from the gluon free energy (3.3) is [48] The hard contribution from the gluon HTL counterterm (3.13) is [48] F (h) The hard contribution from the quark free energy (3.5) is The hard contribution from the quark HTL counterterm (3.15) is Since scalars are bosons, the sum-integrals in (3.6) are the same those used for gluons. After expansion, we obtain the hard contribution to the LO scalar free energy Using the results for sum-integrals contained in the Appendixes B and C of Refs. [48,49], Eq. (5.5) reduces to The scalar HTL counterterm is given in (3.14), after expansion, we get which can be reduced to Note that the first terms in (5.1), (5.3) and (5.6) cancel the order-0 term in the coefficient of mass squared in (5.2), (5.4) and (5.8), respectively.

Soft contributions
The soft contributions to the thermodynamical potential come from the n = 0 Matsubara mode (P 0 = 0) in the resulting bosonic sum-integrals. For fermions, since P 0 = (2n+1)πT = 0 for integer n, the quark momentum is always hard; therefore, quarks do not have a soft contribution. For gluons, in the soft limit P → (0, p), the HTL gluon self-energy functions reduce to Π T (P ) = 0 and Π L (P ) = m 2 D . For scalars, in this limit the propagator reduces to ∆ s (P ) = −1/(p 2 + M 2 D ) where, here p 2 = p 2 . The soft contribution to the gluon free energy (3.3) is The soft contribution from the gluon HTL counterterm (3.13) is The soft contribution from scalar free energy (3.6) is which can be reduced to The soft contribution from the scalar HTL counterterm (3.14) is then it can be reduced as (5.14)

Two-loop sum-integrals
Since the two-loop sum-integrals have an explicit factor of λ, we only need to expand these sum-integrals to order m 3 D m D /T 3 to include all terms through λ 5/2 . Since these integrals involve two momentum integrations we will expand contributions from hard loop momentum and soft loop momentum for each momentum integral. For bosons, this gives three contributions which we will denote as (hh), (hs) and (ss). For fermions, since their momentum is always hard, there will be only two regions (hh) and (hs). In the (hh) region, all three momentum are hard p ∼ T , while in the (ss) region, all the three momentum are soft, p ∼ gT . In the (hs) region, two of the three momenta are hard and the other is soft.

Contributions from the (hh) region
In the (hh) region, the self energies are suppressed by m 2 D /T 2 , M 2 D /T 2 and m 2 q /T 2 , so we can expand in powers of Π T , Π L , Σ, and M 2 D . The (hh) contribution from gluon self energy (3.8) is [ The (hh) contribution from (3.9) can be expanded as where P + R − Q = 0. Using the sum-integral formulas in Appendix C of Ref. [49] this reduces to The (hh) contribution to (3.11) can be expanded as where P + R − Q = 0. Using the sum-integral formulas in Appendix C of Ref. [49] this reduces to (5.20)

Contributions from the (hs) region
In the (hs) region, the soft momentum can be any bosonic momentum. The functions that multiply the soft propagators ∆ T (0, p), ∆ X (0, p), or ∆ s (0, p) can be expanded in powers of the soft momentum p. In terms involving ∆ T (0, p), the resulting integrals over p have no scale and vanish in dimensional regularization. The integration measure p scales like m 3 D for gluon momentum and M 3 D for scalar momentum, respectively. The soft propagators ∆ X (0, p) and ∆ s (0, p) scale like 1/m 2 D and 1/M 2 D , respectively, and every power of p in the numerator scales like m D for gluon momentum and M D for scalar momentum.
The (hs) contribution from the gluonic free energy graphs (3.8) is [48] F (hs) The (hs) contribution from quark self energy (3.10) is Note that the sign on the second term differs from Ref. [49]. This is due to an incorrect sign in the HTL-corrected quark-gluon three vertex in Ref. [49], which we discuss in the Appendix A.9.
For the (hs) contributions to (3.9) and (3.11), like QCD, after expansion there will be terms of contributing at O(g) and higher. For terms that are already of order g 3 , we can set R = Q for soft momentum P . For terms that are O(g), we must expand the sum-integral to second order in p, and then perform the angular integration for p, where the linear terms in p vanish and quadratic terms of the form p i p j can be replaced by p 2 δ ij /d. Therefore, the (hs) contribution from (3.9) can be written as where P + R − Q = 0. Using the sum-integral formulas from Appendix C of Ref. [49] this reduces to The (hs) contribution from (3.11) can be written as where P + R − Q = 0. Again using the sum-integral formulas from Appendix C of Ref. [49] this reduces to F (hs) (5.26)

Contributions from the (ss) region
In the (ss) region, all bosonic momentum are soft, and the gluonic HTL correction functions T P , T 000 , and T 0000 vanish. The gluonic self-energy functions at zero-frequency are Π T (0, p) = 0 and Π L (0, p) = m 2 D . The scales in the integrals come from the gluonic longitudinal propagator ∆ L (0, p) = 1/(p 2 + m 2 D ) and scalar propagator ∆ s (0, p) = −1/(p 2 + M 2 D ). Therefore for bosons, in dimensional regularization, at least one such propagator is required in order for the integral to be nonzero, and there is no (ss) contributions coming from fermionic diagrams.
The (hs) contribution to the gluonic free energy graphs (3.8) is The (ss) contribution to (3.9) can be expanded as , (5.28) where P + R − Q = 0. Again using the sum-integral formulas from Appendix C of Ref. [49] this reduces to

HTL thermodynamic potential
In this section, we calculate the thermodynamic potential Ω(T, λ, m D , M D , m q , δ = 1) explicity, first to LO in the δ expansion and then to NLO.

Leading order
As we mentioned in Sec. 3, the leading order thermodynamic potential is the sum of the contributions from one-loop diagrams and the leading order vacuum energy counterterm. The contributions come from the one-loop diagrams is the sum of (5.1), (5.3), (5.6), (5.9) and (5.12). After multiplying by the appropriate coefficients in (3.2), one obtains where F ideal is the free energy density of N = 4 SYM in the ideal gas limit andm D ,M D , m q andμ are dimensionless variables, defined aŝ Since the leading order vacuum energy counterterm ∆ 0 E 0 should cancel the divergences in the one-loop free energy, we obtain After adding the leading order vacuum renormalization counterterm, our final result for the renormalized LO HTLpt thermodynamic potential is

Next-to-leading order
The next-to-leading order corrections to the thermodynamic potential include all of the two-loop free energy diagrams, the gluon, quark, scalar counterterms in Fig. 3, and the renormalization counterterms. The contributions from the two-loop diagrams include all terms through order g 5 is the sum of (5.  The total NLO HTL counterterm contribution is the sum of (5.2), (5.4), (5.8), (5.10) and (5.14) multiplied by the Casimirs in (3.12) The ultraviolet divergences in (6.5) and (6.6) will be removed by the renormalization of the vacuum energy density E 0 and the HTL mass parameters m D , M D , and m q . The renormalization counterterm contribution at linear order in δ is denoted by ∆Ω NLO in Eqs. (3.1) and (3.16). We cannot obtain its form directly from Ref. [49] due to the fact that there are contributions coming from the scalar fields in N = 4 SYM theory, but we can use the same method as in QCD. The form of ∆ 1 E 0 can be obtained by expanding (2.10) to first order in δ, which is This renormalization counterterm cancels the divergences in Ω HTL (6.6). For the two-loop self energy, as we can see, the divergent terms are q which is a big difference from QCD, we cannot use the following formula This is because we cannot get the two mixed term. In order to cancel the ultraviolet divergence for two-loop self energy, the simplest form is in Eq. (3.16). Using (2.11), (6.7), and (3.16), one finds Adding the leading order thermodynamic potential in (6.4), the two-loop free energy in (6.5), the HTL gluon and quark counterterms in (6.6), and the HTL vacuum and mass renormalizations in (6.10), our final expression for the NLO HTLpt thermodynamic poten-tial in N = 4 SYM is Note that this result reproduces the perturbative expansion given in (1.2) through O(λ 3/2 ) in the weak-coupling limit. This can be verified by takingm D ,M D , andm q to be given by their leading-order expressions (6.16) and truncating the resulting expansion in the 't Hooft coupling at O(λ 3/2 ).

Gap equations
The gluon, scalar, and quark mass parameters m D , M D , and m q are determined by using the variational method, requiring that the derivative of Ω NLO with respect to each parameter is zero The first equation giveŝ The second equation giveŝ  The third equation giveŝ −3M 2 D 4 log 2 + 3γ + 3 logμ 2 + 8m 2 q γ + 2 log 2 + logμ 2 Note that the terms proportional tom 2 q in Eqs. (6.14) and (6.15) can be written in terms ofM D andm D by using (6.13).
In practice, one must solve these three equations simultaneously in order to obtain the gap equation solutions form 2 q (λ),m 2 D (λ), andM 2 D (λ). In Fig. 4 we present our numerical solutions to these three gap equations scaled by the corresponding leading-order weak- In all three panels, the black line is the solution when taking the renormalization scalê µ = 1, the red dashed line isμ = 1/2, and the blue long-dashed line isμ = 2. As can be seen from Fig. 4, the gap equation solution form q does not approach its perturbative limit when λ is approaches zero. This is similar to what was found in NLO HTLpt applied to QCD [49].

Thermodynamic functions
The NLO HTLpt approximation to the free energy is obtained by evaluating the NLO HTLpt thermodynamic potential (6.11) at the solution of the gap equations (6.12) The pressure, entropy density, and energy density can then be obtained using Note that due the conformality of the SYM theory, in all three of these functions, the only dependence on T is contained in the overall factor of F ideal . As a result, when scaled by their ideal limits, the ratios of all of these quantities are the same, i.e. P/P ideal = S/S ideal = E/E ideal .

Numerical results
In  and strong coupling limits [5]. Finally, the grey dotted lines indicate the strong and weak coupling limits of 3/4 and 1, respectively.
As can be seen from Fig. 5, the LO and NLO HTLpt predictions are close to one another out to λ 2. Computing the ratio of the NLO and LO results, we find that they are within ∼ 5% of one another in this range. This is a much smaller change from LO to NLO than is found using the naive weak-coupling expansion. We also observe that the size of the scale variation (shown as shaded red and blue bands) decreases as one goes from LO to NLO. Comparing the bands at λ = 1 we find that the LO HTLpt variation aroundμ = 1 is on the order of 2%, whereas the NLO order HTLpt variation is 0.3%. For λ 6 the NLO HTLpt is below the value expected in the strong coupling limit. At smaller couplings, λ 1 we observe that the NLO HTLpt result is very close to the R [4,4] Padé approximant. This could be coincidental, however, it is suggestive that somehow the R [4,4] Padé approximant may provide a reasonable approximation to N = 4 SYM thermodynamics despite its ad hoc construction.
A similar conclusion was obtained in a prior study of HTL resummation in N = 4 SYM thermodynamics [5]. In their work, Blaizot, Iancu, Kraemmer, and Rebhan (BIKR) used the Φ-derivable framework to obtain an approximately self-consistent approximation to the scaled entropy density. In their approach, the one-loop Φ-derivable result for the entropy density has approximate next-to-leading-order accuracy (NLA), so it should be comparable to the our NLO HTLpt result. In Fig. 6 we present a comparison of our NLO HTLpt result with the BIKR NLA result [5]. In this Figure, the blue line with a blue shaded band is our NLO HTLpt result and the solid red line with a red shaded band is the NLA result from Ref. [5]. The dashed and dotted lines are the same as the previous figure. We find that, forμ = 1, the two calculations are within ≤ 2% of one another for λ 6. We observe that the NLO HTLpt result has a smaller scale variation than the NLA result at all couplings  shown. Finally, in Fig. 7 we present a comparison of all results for λ ≤ 2. The various weakcoupling lines in this Figure are the same as in Fig. 5. As can be seen from this Figure, there is excellent agreement between the NLA calculation of BIKR and NLO HTLpt in this coupling range. We also see that the R [4,4] Padé approximant overlaps with both calculations at smaller λ. Given the agreement between our NLO results and the BIKR NLA results in this range of 't Hooft coupling, one can try to estimate the range of temperatures this might map to in a real-world QGP. This is a fraught endeavor, however, since one can choose to match a variety of quantities with, to our knowledge, no unique prescription. In Ref. [5] the authors advocated matching the scaled entropy density. For the purposes of a ball-park estimate, we will follow their suggestion. State-of-the-art lattice data for the scaled entropy density indicates that corrections to the ideal limit saturate above approximately T ∼ 3T c ∼ 450 MeV at value of S QCD /S QCD,0 ∼ 0.8 − 0.85 [55]. Matching this to the same ratio in N = 4 SYM, one finds from Fig. 6 that this requires λ ∼ 3 − 4 usingμ = 1. 3 Our results suggest that the NLO HTLpt result for N = 4 SYM can be trusted to with high accuracy for λ 2. This provides motivation for extending our calculation to NNLO.

Conclusions
In this paper we have extended the LO and NLO HTLpt calculation of the thermodynamic potential in QCD to N = 4 SYM theory. We have presented results for the LO and NLO HTLpt predictions for the thermodynamics of N = 4 SYM for arbitrary N c . We found that it is possible to extend the range of applicability of perturbative calculations of thermodynamics in N = 4 SYM theory to intermediate couplings, albeit using involved resummations. We compared our NLO HTLpt results to approximately self-consistent resummations obtained previously in Ref. [5] and found them to be in excellent agreement with our NLO HTLpt results for the scaled entropy density for λ 6. Compared to the method used in Ref. [5], our HTLpt results are manifestly gauge-invariant and the HTLpt framework allows for systematic extension of the calculation to higher loop orders. It would be interesting to extend the HTLpt results obtained here to NNLO as has been done in QCD. For this purpose, it seems necessary to first establish the naive perturbative corrections to SUSY thermodynamics at orders g 4 and g 5 . Work along these lines is in progress.

A.1 Gluon polarization tensor
In N = 4 SYM, there are six diagrams that contribute to the LO gluon self energy. In the HTL limit, for massless bosons and fermions, the gluon polarization tensor was derived in Ref. [56] Π µν ab (p) = −g 2 N c δ ab where f (k) ≡ 2n g (k) + 8n q (k) + 6n s (k) is the effective one-particle distribution function for an N = 4 SYM theory. The coefficients of n g , n q , n s are equal to the number of degrees of freedom of the gauge field, fermions, and scalars. Since Π µν ab (p) is symmetric and transverse in its Lorentz indices, it is gauge independent.
We can define the Debye mass for the gauge field using which has been given previously in Refs. [3] and [56]. Then using integration by parts from Ref. [38] applied to (A.1), we obtain where y µ ≡ k µ /|k| = (1, k/|k|) ≡ (1,ŷ). After integration over the length of momentum |k|, the HTL gluon polarization tensor can be written as Where we have introduced a rank-two tensor T µν (p, q) which is defined only when p + q = 0 as T µν (p, −p) = y µ y ν p · n p · y ŷ . (A.5) The angular brackets indicate averaging over the spatial direction of the light-like vector y. The tensor T µν is symmetric in µ and ν, and satisfies the "Ward identity" As a result, the polarization tensor Π µν is also symmetric in µ and ν and satisfies The gluon polarization tensor can also be expressed in terms of two scalar functions, the transverse and longitudinal polarization functions Π T and Π L , defined by

A.5 HTL quark counterterm
The insertion of an HTL quark counterterm into a quark propagator is where Σ(p) is the HTL quark self energy given in (A.27).

A.6 Scalar self-energy
There are four diagrams that contribute to the scalar self energy in N = 4 SYM theory. In the HTL limit, the scalar self energy P AB ab was computed in Ref. [56] for massless bosons and fermions P AB ab (p) = g 2 N c δ ab δ AB After integration over the length of the three-momentum |k|, the HTL scalar self energy reduces to P AA aa (p) = g 2 N c T 2 = λT 2 = M 2 D , (A. 36) where M 2 D is the adjoint scalar mass, which has been given in Ref. [3,56]. We can see that it satisfies M 2 D = m 2 D /2 = 2m 2 q . Note that the scalar self energy is a constant, which means that it only affects the scalar propagator and not the scalar-gluon and scalar-quark vertices in N = 4 SYM theory. This is due to the fact that the HTL Lagrangian density L HTL is a combination of the fields and their corresponding covariantized self energies and the HTL vertices are obtained by expanding the covariant derivatives D µ appearing in the HTL effective Lagrangian in powers of the gauge field A µ . Since there are no covariant derivatives appearing in the scalar contribution to the HTL effective action (2.9), the scalar-gluon vertices will not receive corrections in HTLpt. 4,5

A.7 Scalar propagator
The Feynman rule for the scalar propagator is iδ ab δ AB ∆ s (p) , (A.37) 4 Ref. [38] details the steps necessary to obtain the QCD HTLpt propagators and vertices from the QCD HTL effective action for both equilibrium and non-equilibrium systems. 5 We are grateful for the authors of Ref. [56] for bringing this to our attention.

A.8 HTL scalar counterterm
The insertion of an HTL scalar counterterm into a scalar propagator is −iδ ab δ AB P AA aa (p) .
(A. 40) where P AA aa (p) is the HTL scalar self energy given in (A.36).

A.9 Quark-gluon vertex
The quark-gluon vertex with incoming gluon momentum p, incoming quark momentum r, and outgoing quark momentum q, Lorentz index µ, and color indices a, b, c is Γ µ,ij abc (p, q, r) = −gf abc δ ij γ µ + m 2 qT µ (p, q, r) = −gf abc δ ij Γ µ (p, q, r) . (A.41) Note that the sign on the second term differs from Ref. [49]. This appears to be a typo in the original reference. The rank-one tensorT µ in the HTL correction term is only defined for p + r − q = 0T µ (p, q, r) = y µ y (y · r)(y · q) Note that the overall sign here differs from Ref. [49]. This appears to be a typo in the original reference. The quark-gluon vertex therefore satisfies the Ward identity p µ Γ µ (p, q, r) = S −1 (q) − S −1 (r) . (A.44)

A.13 Scalar-gluon four vertex
The scalar-gluon four vertex is independent on the direction of the momentum, and it can be expressed as Γ µν,AB abcde (p, q, r, s) = −2ig 2 g µν δ AB f ade f bce , (A. 53) and satisfies δ ac δ bd δ AB Γ µν,AB abcde (p, q, r, s) = (2ig 2 )(6N c d A )g µν . (A.54) A.14 Quark-scalar vertex Since fermions have different interactions with the scalar (X p ) and pseudoscalar (Y q ) degrees of freedom, there are two kinds of vertex needed. One is quark-scalar vertex, with incoming scalar momentum p, outgoing quark momentum q and incoming quark momentum r, and their corresponding colors a, b, c respectively. This vertex can be written as