${\cal N}=4$ supersymmetric Yang-Mills thermodynamics to order $\lambda^2$

We calculate the resummed perturbative free energy of ${\cal N}=4$ supersymmetric Yang-Mills in four spacetime dimensions ($\text{SYM}_{4,4}$) through second order in the 't Hooft coupling $\lambda$ at finite temperature and zero chemical potential. Our final result is ultraviolet finite and all infrared divergences generated at three-loop level are canceled by summing over $\text{SYM}_{4,4}$ ring diagrams. Non-analytic terms at ${\cal O}({\lambda}^{3/2}) $ and $ {\cal O}({\lambda}^2 \log\lambda )$ are generated by dressing the $A_0$ and scalar propagators. The gauge-field Debye mass $m_D$ and the scalar thermal mass $M$ are determined from their corresponding finite-temperature self-energies. Based on this, we obtain the three-loop thermodynamic functions of $\text{SYM}_{4,4}$ to ${\cal O}(\lambda^2)$. We compare our final result with prior results obtained in the weak- and strong-coupling limits and construct a generalized Pad\'{e} approximant that interpolates between the weak-coupling result and the large-$N_c$ strong-coupling result. Our results suggest that the ${\cal O}(\lambda^2)$ weak-coupling result for the scaled entropy density is a quantitatively reliable approximation to the scaled entropy density for $0 \leq \lambda \lesssim 2$.


Introduction
The perturbative expansion of the free energy of N = 4 supersymmetric Yang-Mills in four dimensions (SYM 4,4 ) with N c colors and gauge coupling g can be written in the form lim λ→0 F ∼ T 4 a 0 + a 2 λ + a 3 λ 3/2 + a 4 + a 4 log λ λ 2 + O(λ 5/2 ) , (1.1) where λ = g 2 N c is the 't Hooft coupling, which does not run and is independent of the temperature. 1 This expansion is identical in form to the perturbative expansion of the QCD free energy [1][2][3][4][5]. The leading term in this expression is the free energy of an ideal SYM 4,4 plasma and the O(λ) correction can be obtained by computing two-loop Feynman diagrams. Naively, the next contribution is O(λ 2 ) and comes from three-loop contributions. However, a problem emerges because one finds uncanceled infrared divergences at the three-loop level if one uses bare propagators. The same issue occurs in QCD and, in this case, the infrared divergences can be eliminated by summing over the so-called ring diagrams [1][2][3][4]. The solution is the same in SYM 4,4 , with the key difference between QCD and SYM 4,4 being the number and types of degrees of freedom, since the SYM 4,4 theory has six scalar fields and four Majorana fermions. In order to cancel the three-loop infrared divergences, similar to the case of QCD, one can reorganize the perturbative calculation to take into account the soft thermal masses of the gluon and scalar fields [6][7][8]. In QCD, through O(λ 5/2 ), these soft mass scales result in non-analytic contributions to the free energy at orders λ 3/2 , λ 2 log λ, and λ 5/2 [1][2][3][4]9] and an analytic contribution at order λ 2 . In the weak-coupling limit, the free energy of SYM 4,4 has been calculated through order λ 3/2 , with the result being [10][11][12][13] 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. The first three terms in (1.2) map to a 0 = 1, a 2 = − 3 2π 2 and a 3 = 3+ √ 2 π 3 in eq. (1.1). The aim of our work is to compute the coefficients a 4 and a 4 in eq. (1.1).
For this purpose, we will make use of the regularization by dimensional reduction (RDR) method to (1) simplify some aspects of the calculation using dimensional reduction and (2) regulate any intermediate divergent momentum-space integrals encountered during the calculation [14]. The dimensional reduction method was first discussed by Brink, Schwarz, and Scherk for supersymmetric Yang-Mills theories in 2, 4, 6, and 10 dimensions, and was applied to obtain various Yang-Mills theories with extended supersymmetry in 2 and 4 dimensions [15]. Using this method one can show that SYM 4,4 can be obtained by dimensional reduction of 10-dimensional N = 1 supersymmetric Yang-Mills theory (SYM 1,10 ). We will make use of this fact to simplify the calculation of the three-loop SYM 4,4 diagrams by instead computing massless three-loop diagrams in SYM 1,10 followed by dimensional reduction to SYM 4,4 in the RDR scheme. We note that there are inconsistencies in Siegel's RDR scheme which have been pointed out by several authors including Siegel himself [16][17][18][19]. However, these inconsistencies only manifest themselves at high-loop order or under application of Fierz reordering identities and can be pushed to higher orders by making use of the superfield formalism [17]. As discussed in Ref. [17], the breaking of supersymmetry in Siegel's RDR scheme first appears at 3-loop order in the SYM 1,10 beta function and at higher-loop orders for other quantities such as the propagator of the vector supermultiplet. Since a three-loop graph contributing to the running coupling maps to a four-loop vacuum contribution, the three-loop calculation of the thermodynamic potential presented herein should be free from such ambiguities.
In addition to the calculational simplification provided by dimensional reduction, the RDR scheme also provides a regularization method that manifestly preserves supersymmetry. This method was first introduced by Siegel [14] and implements a modified form of dimensional regularization [20][21][22] which manifestly preserves gauge invariance, unitarity, and supersymmetry by making use of dimensional reduction [15,23]. The application of the RDR method to N = 1, 2, and 4 SYM was considered by Avdeev and Vladimirov [17], and the authors presented two equivalent methods for the evaluation of SYM Feynman diagrams based on dimensional reduction. Since Siegel's original paper, other authors have considered variations of RDR in which the size of the representations can depend on the number of spatial dimensions, however, in these variants one must introduce additional degrees of freedom (so-called -scalars) in order to preserve supersymmetry [24]. This gives equivalent results to the original RDR scheme, but complicates the calculation since additional degrees of freedom and Feynman rules must be implemented. For this reason, we make use of Siegel's original RDR scheme. Put simply, in Siegel's RDR scheme, one maintains supersymmetry by keeping the size of the bosonic, fermionic, and scalar representations fixed when the number of spatial dimensions is changed. As a result of this prescription, the cancellations between the bosonic and fermionic degrees of freedom necessary to maintain supersymmetry are automatically preserved.
In the resulting RDR scheme, when computing vacuum contributions to the free energy, one fixes the dimension of the field representations to be integer valued and takes all momentum to be 4 − 2 dimensional vectors, where is an infinitesimal. Although the use of dimensional reduction results in a dramatic simplification in the computation of vacuum contributions since it is more straightforward to calculate higher-loop diagrams in SYM 1,10 than in SYM 4,4 , one limitation of using the dimensional reduction is that, although there is exact equivalence between the results of vacuum graphs computed using dimensional reduction of 10-dimensional SYM 1 theory and the results of vacuum graphs computed directly in SYM 4,4 , it is not possible to treat soft contributions to the free energy in the same manner. For this purpose, one must introduce dressed gauge and scalar propagators directly in SYM 4,4 in order to resum the supersymmetric ring diagrams and eliminate infrared divergences generated at three-loop order. The dimensional reduction method was first used for the computation of two-loop SYM 4,4 thermodynamics in ref. [12] where it was shown that, with this method, one can compute the hard contributions using dimensionally reduced SYM 1,10 and the soft contributions directly in SYM 4,4 . We will adopt the same strategy: (1) use dimensional reduction applied to SYM 1,10 to compute the massless three-loop vac-uum graphs and (2) compute the necessary soft contributions and counterterms directly in SYM 4,4 by dressing the gauge and scalar propagators.
Finally, we note that in the opposite limit of strong coupling, the behavior of the SYM 4,4 free energy has been computed using the anti-de Sitter space/CFT (AdS/CFT) correspondence [25]. In the large-N c limit one has [26] where, in the large-N c limit, it is expected that only powers of λ −3/2 appear in the strongcoupling expansion. 2 In the results section we use this information to constrain the coefficients appearing in a large-N c generalized Padé approximant constructed from terms of the form λ n/2 and λ n/2 log λ with n being a positive integer. We demonstrate that the resulting generalized Padé is free from singularities at all λ, reproduces the weak-coupling expansion through O(λ 2 , λ 2 log λ), and reproduces the strong-coupling expansion to all known orders, with no terms containing log λ through O(λ −5/2 log λ). The structure of our paper is as follows. We begin with a brief summary of SYM 1,10 and its relation to SYM 4,4 in the RDR scheme. In sec. 3, we recall the basics of SYM 4,4 and introduce the resummations necessary to obtain the three-loop corrections to the free energy. In sec. 4, we list all contributions to the free energy up to three-loop order and describe the key steps necessary to obtain our final results. In this section we also present our final results at each loop order truncated at O( 0 ). In sec. 5 we present our final result for the O(λ 2 ) and O(λ 2 log λ) coefficients, present comparisons with past results, and compare to a generalized Padé approximant constructed using the new constraints provided in the weak-coupling limit. Finally, in sec. 6 we present our conclusions and an outlook for the future.

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.
2 Dimensional reduction of SYM 1,D in the RDR scheme The action of N = 1 supersymmetric Yang-Mills in D dimensions (SYM 1,D ) can be written in Minkowski space as [12,15] where M, N = 0, · · · , D − 1, and the field strength tensor is is the covariant derivative in the adjoint representation of SU (N c ). 3 The definition of the gauge field is the same as QCD and A M can be 2 If one considers sub-leading corrections in 1/Nc in the strong-coupling limit, the entropy density can contain additional fractional powers of λ [27]. It is also possible that there could be logarithms of the 't Hooft coupling induced by massless gravity modes in this case. 3 In practice, the gauge group can be any semi-simple Lie group. Since our target theory is SYM4,4 with the SU (Nc) group, here we also use the SU (Nc) group.
expanded as A M = A a M t a , with real coefficients A a M , 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 metric tensor is g M N = (1, −1, · · · − 1) and the Γ-matrices satisfy with I n being an identity matrix of dimension n. SYM theories with a different number of supercharges # SC can be obtained by choosing the appropriate value of D for which the number of supercharges is maximal In order to maintain supersymmetry, the number of bosonic and fermionic degrees of freedom must be equal, which implies that the fermions in D max = 10 must satisfy the Majorana-Weyl condition, while they satisfy the Weyl condition in D max = 6 and D max = 4 [12]. Because of these constraints, in general, the number of degrees of freedom is equal to D max − 2. In this way, one learns that the n in (2.3) should be the dimension of the spinors of the maximal SYM theory, and equal to D max − 2, which implies Tr I n = D max − 2.
To quantize the theory, gauge-fixing and ghost terms should be added to the Lagrangian density. In general covariant gauge they are the form with ξ being the gauge parameter. In general, we will be interested in SYM theories with # SC supercharges in dimensions D ≤ D max , where D is an integer. Any D ≤ D max dimensional SYM theory can be obtained by dimensional reduction from the corresponding maximal N = 1 SYM theory in D = D max [15]. To connect these different theories in a manifestly supersymmetric manner one can use a scheme called regularization by dimension reduction (RDR) [14,17,24,28]. In this scheme, the D-dimensional space is decomposed into a direct sum of D-and (D − D)dimensional subspaces. One can obtain N = 1 in D = 8, N = 2 in D = 6, N = 4 in D = 4 and N = 8 in D = 2 theories starting from N = 1 with D max = 10. One can obtain N = 2 in D = 4 and N = 4 in D = 2 starting from N = 1 with D max = 6, while one can obtain N = 2 in D = 2 starting from D max = 4 [12,15]. In the following, we will use d to represent the dimension of the momentum in all theories.
The evaluation of Feynman diagrams for theories that are obtained by dimensional reduction of SYM 1,D can be carried out in two equivalent ways, one is working throughout in (D − 2 ) ⊕ (D − D + 2 ) space for the theories we target, with being an infinitesimal quantity used for regularization and taking the dimension of the momentum-space to be d = D−2 . This prescription results in one having to introduce -scalars into the Lagrangian in order to preserve supersymmetry [24]. A simpler way to maintain supersymmetry is to take all fields in (2.1) to be D-dimensional tensors or spinors and all momentum to be d = D − 2 vectors. We will use this second RDR scheme since it is the most transparent and efficient for computations [14].
With this in mind, up to thermal mass corrections, for any SYM theory that is characterized by (D max , d), massless (unresummed) contributions to the free energy can be calculated perturbatively by computing vacuum Feynman diagrams using N = 1 SYM in D = D max and then restricting the momentum in loops to d dimensions to obtain the result in the target theory.  By decomposing the vector field in 10 dimensions, one obtains six independent real scalar fields which are represented by a multiplet where X p and Y q are 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 the vector Φ. Therefore Φ A , X p , and Y q can be expanded as The action and Lagrangian that generates the perturbative expansion for SYM 4,4 in Minkowski-space can be expressed as where µ, ν = 0, 1, 2, 3 and α p and β q are 4 × 4 matrices that satisfy Up to unitary transformation, their explicit form is 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.
Since we want to obtain the thermodynamic functions up to O(λ 2 ), we need to calculate Feynman diagrams through three loop order. However at three loop level in QCD [1,2], infrared divergences appear that need to be canceled by summing over the ring diagrams appearing in the thermal mass counterterm. As a result, there will be non-analytic contributions to the free energy at order λ 3/2 and λ 2 log λ, representing a breakdown of naive perturbation theory due to infrared thermal physics. In QCD, the need for such reorganizations of perturbation theory stems from the behavior of the thermal propagator at soft momentum, p soft ∼ √ λT . In the weak-coupling limit, the temperature is much higher than p soft and one finds that the thermal mass cannot be treated as a perturbation in the static propagator of bosonic fields.
As detailed in refs. [1,2], in order to systematically resum the necessary diagrams, we need to modify the static bosonic propagators by incorporating gluon and scalar thermal masses, m D and M , respectively. Such a reorganization is necessary in order to make finitetemperature perturbation theory well-defined beyond O(λ). In practice, the thermal masses are determined at leading order by computing the gluon and scalar self-energies Π µν and P at zero four-momentum. In this manner, a scalar thermal mass is generated for Φ and a gluonic thermal mass is generated for the A 0 but not for A. We will introduce m D and M only for the zero Matsubara modes of the gluon and scalar fields, generalizing the method introduced by Arnold and Zhai [1,2]. In the context of finite temperature QCD, such resummations have been carried out through O(λ 3 log λ) [1][2][3][4][5] and reorganizations of perturbation theory based on hard thermal loop perturbation theory and the Φ-derivable approach have shown that one can extend the radius of convergence of weak-coupling treatments to intermediate couplings g ∼ 2 [34][35][36][37][38][39][40]. We note that such perturbative reorganizations have also been carried out to two-loop order in SYM 4,4 [41,42]. Following Arnold and Zhai, in this work we introduce thermal masses, m D and M , only for the zero Matsubara modes of the gluon and scalar fields. The resulting reorganized Lagrangian density in frequency space can be rewritten as where δ p 0 is shorthand for the Kronecker delta function δ p 0 ,0 . Then we absorb the two A 2 0 and Φ 2 terms in the curly brackets into our unperturbed Lagrangian L 0 , and treat the two terms outside the curly brackets as a perturbation.
At the n th loop order in the Arnold and Zhai reorganization the leading term contributing to the thermodynamic potential is O(λ n−1 ); however, due to the thermal masses, each loop order can contain terms that are higher order in λ. Such corrections make the map between loop order and the contributing perturbative orders more complicated. Note importantly that it is not possible for the thermal mass corrections to generate lower-order terms, they only generate higher-order terms than O(λ n−1 ) at each loop order. As a result, there can be no modification of the O(λ 2 ) contribution coming from four or higher-order loop graphs. We also note that terms containing logarithms formally contribute at the same order as the power of λ multiplying them, i.e. λ 2 and λ 2 log λ both contribute at O(λ 2 ).
Similar to QCD, in order to obtain the thermodynamics up to O(λ 2 ), we need to calculate massive Feynman diagrams up to two-loop order and we can take the three-loop diagrams to have bare propagators. Since the propagators in the three-loop diagrams are massless, one can use the RDR method to compute these contributions in SYM 1,10 and then using g M M = 10, D max = 10, and Tr I n = 8, with all momentum integrals computed in d = 4 − 2 . This dramatically reduces the number of diagrams that one must evaluate. For the massive two-loop diagrams, however, we must compute directly in SYM 4,4 , since we must dress the gluon and scalar propagators differently.
As explained in ref. [1,2], the one-and two-loop free energies are gauge parameter independent in the resummed scheme in QCD. In SYM 4,4 , one has additional scalars and fermions, however, their propagators do not depend on the gauge parameter. Since the SYM 4,4 gluon propagator has the same form as in QCD (up to the definition of the gluon thermal mass m D ), the one-loop free energies are still gauge invariant. In addition, we will demonstrate that the two-loop diagrams which contain a three/four-gluon vertex, a gluon-ghost vertex, and/or gluon-quark vertex are gauge invariant in SYM 4,4 .

NNLO thermodynamics of SYM 4,4
Since the RDR method is based on analytically continuing only the number of coordinates and the momentum, the regularized momentum integral can be defined in a similar manner as in QCD. Since we are considering SYM 4,4 , the regularized momentum integral in d = 4 − 2 dimensions can be defined as where withμ being an arbitrary momentum scale. It is also convenient to introduce various invariants associated with the representations of the SU(N c ) gauge group. Denoting the generators of the adjoint representation as (T a ) bc = −if abc and generators of the fundamental representation as t a we define the following group theory factors With the standard normalization

The resummed one-loop free energy
For the one-loop contributions, using the resummed Lagrangian density (3.7), one finds that the one-loop free energy for gluons and scalars will result in a perturbative contribution at O(λ 3/2 ). The one-loop free energy can be written as where d A = N 2 c − 1 is the dimension of the adjoint representation. The Feynman diagrams corresponding to each term are presented in fig. 1. We note that there are four independent Majorana fermions in the adjoint representation giving d F = 4d A . For the scalars, one has d S = 6d A .
Using the resummed gluonic propagator, F 0a can be expressed as where D is the number gluon polarization states. In vacuum, two of the polarizations are unphysical and are canceled by the ghost contribution. We note that this same form can be obtained from the one-loop gluonic hard-thermal-loop perturbation theory (HTLpt) result [42] by setting the gluonic transverse self-energy to zero and the longitudinal selfenergy to m 2 Making use of the resummed scalar propagator, F 0c can be expressed as Since there are no thermal mass contributions for fermions and ghosts, their one-loop free energies have simpler forms with f 0 defined as We note that the resulting one-loop gluon and scalar free energies are equal to the corresponding summation of the leading-order one-loop hard and soft contributions in HTLpt when taking m q → 0 [42]. By combining (4.5), (4.7), (4.8), and (4.9) one obtains (4.12)

The resummed two-loop free energy
The SYM 4,4 two-loop free energy can be written as The Feynman diagrams corresponding to each term are presented in fig. 2. Making use of the resummed gluon and scalar propagators, one obtains The fundamental integrals appearing above are defined as and One finds that the contributions F 1a , F 1b , and F 1c are the same as in QCD [1] and that F 1d is two times the result obtained in QCD in [2], up to the definition of m D and D. As in QCD, the two-loop free energy is gauge independent using our resummation method, so that F 1d and the summation of F 1a , F 1b and F 1c are gauge parameter independent. The only remaining contributions which depend on the gauge parameter are F 1e and F 1g . By using the full propagator in A.1, the gauge terms in F 1e and F 1g are , respectively, and hence cancel. This proves that our resummed two-loop result is gauge-independent in the RDR scheme.
The form of F 1g is similar to F 1b . The first term in each is the result when m D and M are set to zero. The terms containing δ 1 and δ 3 are the corrections in which only one of the three momenta possesses a zero mode (hard-hard-soft). The terms containing δ 2 and δ 4 stem from contributions in which all the three momenta possess zero modes (soft-soft-soft). For F 1g where the quantity in brackets on the second line corresponds to the sum-integral on the first line, but with thermal masses taken to zero. The factors of 2P +Q and Q 0 appearing in the expression come from the gluon-scalar vertex. Making use of the relation (2P + Q) 2 = 2P 2 + 2(P + Q) 2 − Q 2 , δ 3 can be reduced to which can be written in the same form as in eq. (4.20). An expression for δ 4 can be calculated similarly. The form of F 1d and F 1h are similar with the exception of the appearance of I f resum . This is due to the tensor nature of the gluon propagator. The contribution F 1h can be written as . It can be reduced to (4.32) Since the last term in eq. (4.32) is divergence free, one can ignore the M 2 in the denominator. As a result, the above equation can be written in the same form as eq. (4.21). From ref. [1,2] one has and P δ p 0 One finds that the two-loop contributions proportional to λ 3/2 cancel since m 2 D = 24b 1 λ and M 2 = 12b 1 λ. Combining all contributions, we obtain the two-loop resummed free energy (4.37)

The resummed three-loop free energy
As mentioned earlier, the calculation of the massless three-loop vacuum Feynman diagrams in SYM 4,4 can be accomplished more simply in the corresponding SYM 1,10 theory. As a result of this equivalence, one can consider the much smaller set of SYM 1,10 graphs presented in fig. 3, which are topologically equivalent to three-loop QCD vacuum graphs. The threeloop results in SYM 4,4 can be obtained by imposing D = D max = 10, d = 4 − 2 in the SYM 1,D theory.
where in all contributions we have taken d = 4 − 2 . Above 51) and (4.54) The results for I bb ball , I ff ball , I bf ball , and H 3 truncated at O( 0 ) have been given in refs. [1,2]  and The derivations of I bb SYM 1,D , I bf SYM 1,D , and I ff SYM 1,D are presented in app. D. Infrared divergences will be generated by some of the diagrams in fig. 3 due to threemomentum integrations involving massless gluon and scalar propagators in d = 4 − 2 dimensions. These divergences will be canceled by the thermal mass counterterm diagrams in fig. 4. The first diagram in fig. 4 represents the gluonic thermal counterterm contribution, where the shaded blob can be expressed as and Π µν (P ) = Π b µν (P ) + Π f µν (P ), as defined in eqs. (F.1) and (F.2). The form of the second term is equivalent to the form presented in the ref. [1,2], however, we choose the above form for clarity and calculational ease. Similarly, the shaded blob in the second diagram in fig. 4 stands for the scalar thermal counterterm and can be expressed as ∆P(P ) ≡ P(P ) + M 2 δ P 0 . (4.60) The first diagram in fig. 4 is simple and is proportional to (4.64) To simplify further one can use eq. (D.10) We note that this result holds for all N c . As can be seen from this expression, one finds non-vanishing contributions at O(λ 2 ) and O(λ 2 log λ), as anticipated. Equation (5.1) is manifestly finite due to an explicit cancellation between three-loop infrared singularities and the three-loop counterterm diagrams. These cancellations remove all infrared divergent contributions. In addition, there are no remaining poles due to ultraviolet divergences, since the coupling does not run in SYM 4,4 and, hence, no coupling constant renormalization counterterm is required. Based on the result above, we see that the series organizes itself naturally as an expansion in λ/π 2 , suggesting that the weak coupling expansion will break down for λ π 2 . The presence of the logarithm at order λ 2 stems directly from the dressing of the A 0 and scalar propagators. The pressure, entropy density, and energy density can be obtained from (5.1) using standard thermodynamic identities Note that due the conformality of SYM, 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 . In Fig. 5 we plot the scaled entropy density as a function of λ. The green dotted, red dashed, and blue long-dashed curves correspond to the perturbative result truncated at This suggests that adding each perturbative order extends the estimated range of validity in λ by an order of magnitude. In comparison to the convergence of the perturbative QCD free energy we observe that the O(λ 2 ) truncation in SYM 4,4 has P/P ideal = S/S ideal < 1 for λ 10, whereas the O(λ 2 ) truncation in QCD has P > P ideal for λ 3.5 for the central value of the renormalization scale. In contrast, lattice QCD measurements of the pressure find P < P ideal . In QCD, one has to include all contributions through O(λ 5/2 ) in order for the pressure to be less than the ideal pressure at large coupling. This suggests that the perturbative expansion of the SYM 4,4 free energy might have better convergence than the perturbative expansion of the QCD free energy.

Conclusions and outlook
In this paper we computed the thermodynamic function of SYM 4,4 to O(λ 2 ). Our final result, presented in eq. (5.1), extends our knowledge of weak-coupling SYM 4,4 thermodynam-ics to include terms at O(λ 2 ) and O(λ 2 log λ). The appearance the O(λ 2 log λ) contribution can be traced back to non-analytic terms generated due to dressing of the gluon and scalar propagators at finite temperature (ring resummation). All results presented here made use of the RDR scheme which manifestly preserves supersymmetry. Having obtained the O(λ 2 ) and O(λ 2 log λ) coefficients in the SYM 4,4 free energy, we then constructed a large-N c Padé approximant that interpolates between the weak-and strong-coupling limits. The resulting singularity-free Padé approximant (G.1) reproduces the weak-coupling limit through O(λ 2 , λ 2 log λ) and the coefficients and analytic structure of the large-N c strong-coupling limit through O(λ −5/2 ), with no terms containing log λ appearing through this order.
In the near future we plan to also compute the coefficient of λ 5/2 in the SYM 4,4 free energy. For this purpose, one can either adapt the methods presented originally by Zhai and Kastening in QCD [3] to SYM 4,4 or one could consider applying effective field theory methods similar to refs. [4] and [13]. We plan to pursue the first option in the near term due the fact that the results presented herein followed the Arnold, Zhai, and Kastening calculational framework. It would certainly be interesting to consider the application of effective field theory methods to SYM 4,4 thermodynamics since, as emphasized herein, the necessary three-loop diagrams can be evaluated using dimensional reduction of the simpler SYM 1,10 theory. Finally, we mention that we also plan to pursue a three-loop HTLpt calculation of SYM 4,4 thermodynamics in order to further improve the convergence of the successive approximations to the thermodynamic functions in the weak-coupling limit. This would extend our prior two-loop HTLpt calculation of SYM 4,4 thermodynamics to threeloop order [42].
The future work outlined above will extend the perturbative calculation of the SYM 4,4 free energy to the highest order possible before genuinely non-perturbative effects related to magnetic scales in SYM enter [13]. Like QCD, in SYM 4,4 at O(λ 3 ) there will be both perturbative and non-perturbative contributions. The non-perturbative contribution at O(λ 3 ) can be computed on the lattice using effective field theory methods, however, computation of the four-loop perturbative contributions at this order presents a technical challenge.

Note added
In arXiv version 3, we corrected an error in eq. 2) compared to version 2. These changes has been verified by computing the free energy using effective field theory methods. Importantly, our conclusions are unchanged since the correction is numerically small, however, an erratum correcting this error has been submitted to JHEP. We have also taken this opportunity to correct some typos. Based on the Lagrangian density in eq. (3.7) one sees that only the gluon and scalar propagator are modified using this resummation method and, as a consequence, there will be two thermal mass counterterm vertices which must be taken into account.

A.1 The resummed gluon propagator
The Feynman rule for the resummed gluon propagator is where the tensor-valued gluon propagator ∆ µν depends on the choice of gauge fixing. In covariant gauge it can be expressed as

A.2 The resummed scalar propagator
The Feynman rule for the resummed scalar propagator is Note that this is the same form as the λφ 4 resummed propagator presented in ref. [8].

A.3 The counterterm vertex
The gluonic counterterm vertex is The scalar counterterm vertex is and we have taken D = D max . One can obtain (Π ab M N (P )) SYM 1,D directly from the corresponding result in QCD. From the Lagrangian, in the bosonic sector the only difference between SYM 1,D and QCD is the dimension of the metric tensor and momenta. As a result, one can immediately generalize the one-loop self-energy results from QCD. In the fermionic contributions an additional difference stems from the different fermionic spinor and color representations in the two theories, which causes Tr I n = D − 2 and S F = N c . 4 In the hard-thermal-loop limit (HTL), by using the contour integral method and integration by parts, one obtains the bosonic thermal mass in SYM 1,D where d = D. We note that eq. (C.5) is the same as the result that would be obtained by using −δ M N Π M N , but for ease of application of the RDR scheme, we choose the form in eq. (C.5).
In the same way, by calculating the one-loop fermionic self-energy in the HTL limit, one obtains the quark thermal mass By imposing D max = 10 and d = 4 − 2 in (C.5) and (C.6), one obtains Note that these are the same as the bosonic and fermionic thermal masses with D = 4 in SYM 4,4 . However, in the SYM 1,10 theory one cannot obtain the thermal mass for scalars straightforwardly since they do not appear as fundamental degrees of freedom. As a consequence, one cannot use the RDR method based on massive SYM 1,10 diagrams and one must perform the necessary resummations directly in SYM 4,4 .
where D = D max . Inserting (C.1) and (C.2) and expanding, one obtains three terms we define asĪ bb where H 4 , H 5 and H 6 are defined by (D.9) The final results for these integrals are given in app. E. Simplifying and using is to use the method introduced in ref. [2], which is to use the final results for I bb sqed , I bf sqed , and I ff sqed , defined by where ∆Π µν (P ) ≡Π µν (P ) −Π ρρ (0)δ µ0 δ ν0 δ p 0 . For the calculation of H 4 , using eq. (F.4) in eq. (E.3) and simplifying, one obtains From eq. (E.4) one finds I ff sqed is independent of the thermal mass contribution because I f resum = 0. This will not occur in I bb sqed and I bf sqed since I b resum = 0. By using eqs. (E.1), (4.50), (E.2), and (4.51), the final results for H 5 and H 6 can be obtained The one-loop gauge bosonic self-energy of (Π ab µν (P )) SYM 4,4 shown in Fig. 7 contains two parts, Π f,ab µν (P ) = 2λδ ab 2Π f µν (P ) − 2(P 2 δ µν − P µ P ν )

G Large-N c generalized Padé approximant
With new perturbative coefficients in hand one can produce an updated Padé approximant along the lines presented in refs. [11,41]. For this purpose one can make use of the fact that, in the large-N c limit, the strong-coupling expansion becomes a series in powers of λ −3/2 only. Using this information, one can construct a constrained large-N c interpolating function that connects the weak coupling expansion and the strong-coupling expansion. We note, however, that if one considers sub-leading corrections in 1/N c in the strong-coupling limit, the entropy density can contain additional fractional powers of λ [27] and, in this case, there may also be logarithms induced by massless gravity modes. Since knowledge of such corrections is limited, we focus on constructing a Padé approximant that is valid in the large-N c limit only. 5 In order for a Padé approximant to approach a constant in the strong-coupling limit it must be a symmetric or balanced Padé approximant, in which the maximal powers of λ appearing in the numerator and denominator are the same. Based on the large-N c structure of the strong-coupling expansion, we find that the following form can reconstruct all known coefficients in both the weak-and strong-coupling limits S S ideal = 1 + aλ 1/2 + bλ + cλ 3/2 + dλ 2 + eλ 5/2 1 + aλ 1/2 +bλ + 4 3 cλ 3/2 + 4 3 dλ 2 + 4 3 eλ 5/2 , (G.1) 5 Exact solutions for the entropy density in 2+1D O(N ) conformal field theories in the large-N limit have a similar analytic structure, namely logarithms of the coupling appearing in the weak-coupling limit, but not in the strong-coupling limit [43]. We note that, at the order used, namely through O(λ 5/2 ) in the numerator and denominator, the coefficients in the Padé approximant (G.1) obtained are unique. That said, it is possible to add additional terms in both the numerator and denominator, e.g. extending to including O(λ 3 ) terms, which would introduce new coefficients that are unconstrained due to limited information from the weak-and strong-coupling expansions. Based on current information, our final result (G.1) represents the highest order Padé approximant for which all coefficients can be uniquely constrained.