Three-loop HTLpt thermodynamics at finite temperature and chemical potential

We calculate the three-loop thermodynamic potential of QCD at finite temperature and chemical potential(s) using the hard-thermal-loop perturbation theory (HTLpt) reorganization of finite temperature and density QCD. The resulting analytic thermodynamic potential allows us to compute the pressure, energy density, and entropy density of the quark-gluon plasma. Using these we calculate the trace anomaly, speed of sound, and second-, fourth-, and sixth-order quark number susceptibilities. For all observables considered we find good agreement between our three-loop HTLpt calculations and available lattice data for temperatures above approximately 300 MeV.


Introduction
Quantum chromodynamics (QCD) describes the propagation and interaction of quarks and gluons which are believed to be the fundamental constituents of all hadronic matter. Based solely on the QCD Lagrangian it is possible to calculate the finite temperature and chemical potential partition function of QCD which results in the so-called equation of state (EoS). The determination of the QCD EoS is extremely important to the phenomenology of the quark-gluon plasma (QGP). At this time, the most reliable method to calculate the QCD thermodynamic functions at finite temperature and zero chemical potential is lattice gauge theory (see e.g. [1][2][3][4][5][6][7][8][9][10][11][12][13][14]). Importantly, lattice QCD can be used to probe the behavior of QCD matter near the transition temperature where QCD matter undergoes a phase transition from the hadronic phase to the deconfined QGP phase. Near the phase transition, the running coupling is large and non-perturbative methods like lattice QCD must be used. Finite temperature lattice QCD calculations are now quite sound; however, due to the sign problem, it is not straightforward to extend such calculations to finite baryon chemical potential. In practice, it is possible to obtain information about the behavior of the thermodynamic functions at small baryon chemical potential by making a Taylor expansion of the partition function around µ B = 0 and extrapolating the result. This requires the calculation of various quark-number susceptibilities evaluated at zero chemical potential.
Since extrapolations based on a finite number of Taylor coefficients can only be trusted within the radius of convergence of the expansion, it would be nice to have an alternative framework for calculating the finite temperature and chemical potential QCD thermodynamic potential and associated quantities. This is important in light of the ongoing beam energy scan at the Relativistic Heavy Ion Collider (RHIC) and the forthcoming experiments at the Facility for Antiproton and Ion Research (FAIR). As an alternative to lattice QCD calculations, one natural option is to compute the thermodynamic potential using perturbation theory. In principle, this should work since, at sufficiently high temperature, the value of the strong coupling constant is small; however, one does not know a priori how large the temperature should be for this method to result in a good approximation to reality. The calculation of the thermodynamic potential using the weak-coupling expansion in the strong coupling constant, g, has a long history [15][16][17][18][19][20][21][22][23] and the perturbative expansion of the pressure of QCD at both zero [24] and non-zero chemical potential [25][26][27] are now known through order g 6 ln g.
In this paper we calculate the thermodynamic potential at finite temperature and chemical potential to three-loop order in HTLpt. The three-loop thermodynamic potential is renormalized using only known vacuum, mass, and coupling constant counterterms and the final result is completely analytic. The resulting analytic thermodynamic potential is then used to obtain expressions for the pressure, energy density, entropy density, trace anomaly, speed of sound, and various quark number susceptibilities. We find that there is good agreement between our NNLO HTLpt results and lattice data down to temperatures on the order of 300 MeV.
The paper is organized as follows. In section 2 we specify the HTLpt calculational framework and the necessary counterterms to renormalize HTLpt. In section 3 we discuss the diagrams that contribute to the HTLpt thermodynamic potential through NNLO. In section 5 we discuss the mass prescription for the in-medium masses m D and m q . We present our results for the thermodynamic functions and compare them with results from lattice gauge simulations in section 6. In section 7 we present our results for the second-, fourth-, and sixth-order baryon and quark number susceptibilities. We also compare our results for these quantities with available lattice data. In section 8 we summarize and conclude. In appendix A the necessary diagrams are reduced to scalar sum-integrals and expanded in powers of m D /T and m q /T . We list the necessary non-trivial sum-integrals and integrals in appendices B and C. Finally, in appendix D we list some properties of the ℵ functions which appear repeatedly in finite density calculations.

Hard-thermal-loop perturbation theory
The QCD Lagrangian density in Minkowski space can be written as The term ∆L QCD contains the counterterms necessary to cancel ultraviolet divergences in perturbative calculations. The ghost term L gh depends on the form of the gauge-fixing term L gf . In this paper we work in general covariant gauge where L gf = −ξ −1 Tr (∂ µ A µ ) 2 with ξ being the gauge-fixing parameter.
Hard-thermal-loop perturbation theory is a reorganization of in-medium perturbation theory for QCD. The HTLpt Lagrangian density can be written as where the HTL improvement term is [89] where y µ = (1,ŷ) is a light-like four-vector withŷ being a three-dimensional unit vector and the angular bracket indicates an average over the direction ofŷ. The two parameters m D and m q can be identified with the Debye screening mass and the thermal quark mass, respectively, and account for screening effects. HTLpt is defined by treating δ as a formal expansion parameter. By coupling the HTL improvement term (2.3) to the QCD Lagrangian (2.1), HTLpt systematically shifts the perturbative expansion from being around an ideal gas of massless particles to being around a gas of massive quasiparticles which are the appropriate physical degrees of freedom at high temperature and/or chemical potential.
The HTLpt Lagrangian (2.2) reduces to the QCD Lagrangian (2.1) if we set δ = 1. Physical observables are calculated in HTLpt by expanding in powers of δ, truncating at some specified order, and then setting δ = 1. This defines a reorganization of the perturbative series in which the effects of m 2 D and m 2 q terms in (2.3) are included to leading order but then systematically subtracted out at higher orders in perturbation theory by the δm 2 D and δm 2 q terms in (2.3). To obtain leading order (LO), next-to-leading order (NLO), and next-to-next-leading order (NNLO) results, one expands to orders δ 0 , δ 1 , δ 2 , respectively. Note that HTLpt is gauge invariant order-by-order in the δ expansion and, consequently, the results obtained are independent of the gauge-fixing parameter ξ.
If the expansion in δ could be calculated to all orders, the final result would not depend on m D and m q when we set δ = 1. However, any truncation of the expansion in δ produces results that depend on m D and m q . As a consequence, a prescription is required to determine m D and m q as a function of T , µ and α s . Several prescriptions had been discussed in [52] at zero chemical potential. The HTLpt expansion generates additional ultraviolet divergences. In QCD perturbation theory, renormalizability constrains the ultraviolet divergences to have a form that can be cancelled by the counterterm Lagrangian ∆L QCD . We will demonstrate that the renormalization of HTLpt can be implemented by including a counterterm Lagrangian ∆L HTL among the interaction terms in (2.3). There is no all-order proof that the HTL perturbation expansion is renormalizable, so the general structure of the ultraviolet divergences is unknown. However, as shown previously in refs. [44][45][46]52], it is possible to renormalize the NNLO HTLpt thermodynamic potential using only a vacuum counterterm, a Debye mass counterterm, a fermion mass counterterm, and a coupling constant counterterm. The necessary counterterms for renormalization of the NNLO thermodynamic potential are where, with the standard normalization, the QCD Casimir numbers are Note that the coupling constant counterterm (2.7) is consistent with one-loop running of α s .
In practice, in addition to the δ expansion, it is also necessary to make a Taylor expansion in the mass parameters scaled by the temperature, m D /T and m q /T , in order to obtain analytically tractable sum-integrals. An added benefit of this procedure is that the final result obtained at NNLO is completely analytic. In order to truncate the series in m D /T and m q /T one treats these quantities as being O(g) at leading order, keeping all terms that naively contribute to the thermodynamic potential through O(g 5 ). In practice, such an truncated expansion works well [32,90] and the radius of convergence of the scaled mass expansion seems to be quite large, giving us confidence in this approximate treatment of the necessary sum-integrals.

Contributions to the HTLpt thermodynamic potential through NNLO
The diagrams needed for the computation of the HTLpt thermodynamic potential through NNLO can be found in figures 2 and 3 of ref. [52]. In ref. [52] the authors computed the NNLO thermodynamic potential at zero chemical potential. Here we extend the NNLO calculation to finite chemical potential. 1 For this purpose, one needs to only consider diagrams which contain at least one quark propagator; however, for completeness we also list the purely gluonic contributions below. In the results we will express thermodynamic quantities in terms of two dimensionless variables:m D = m D /(2πT ) andμ = µ/(2πT ).
The complete NNLO HTLpt thermodynamic potential can be expressed in terms of these diagrams as where the necessary counterterms at any order in δ can be calculated using eqs. (2.4)-(2.7).
The expressions for the one-and two-loop diagrams above can be found in refs. [44,45]. The expressions for the three-loop bosonic diagrams F g 3a -F g 3m are presented in section 3 of ref. [49], and the three-loop diagrams with fermions F f 3a -F f 3i can be found in section 3 of ref. [50]. The three-loop diagrams specific to QCD, i.e. the non-Abelian diagrams involving quarks, are given by 3) Tr Γ αβ (P, −P, Q, Q)S(Q) ∆ αµ (P )∆ βν (P )Π µν g (P ) ,

NNLO HTLpt thermodynamic potential
One can evaluate the sum-integrals necessary analytically by expanding in the ratios m D /T and m q /T . For details concerning this expansion and intermediate results, we refer the reader to appendix A. We consider first the case that all quarks have the same chemical potential After presenting the steps needed for this case, we give the general result with separate chemical potentials for each quark flavor.

NNLO result for equal chemical potentials
When all quarks have the same chemical potential µ i = µ = µ B /3 we can straightforwardly combine the results for the various sum-integrals. In this case, the unrenormalized threeloop HTLpt thermodynamic potential is where Ω 0 = −d A π 2 T 4 /45. The sum of all counterterms through order δ 2 is where ∆Ω YM is the pure-glue three-loop HTLpt counterterm [52] ∆Ω Adding the total three-loop HTLpt counterterm (4.3) to the unrenormalized three-loop HTLpt thermodynamic potential (4.2) we obtain our final result for the NNLO HTLpt thermodynamic potential in the case that all quarks have the same chemical potential where Ω YM NNLO is the NNLO pure glue thermodynamic potential [49] 2 Note that the full thermodynamic potential (4.5) reduces to thermodynamic potential of ref. [52] in the limit µ → 0. In addition, the above thermodynamic potential produces the correct O(g 5 ) perturbative result when expanded in a strict power series in g [25,26]. 3 2 Note that chemical potential dependence also appears in pure glue diagrams from the internal quark loop in effective gluon propagators and effective vertices. This chemical potential dependence enters through the chemical potential dependence of the Debye mass. 3 There is a mismatch in one term proportional to s2F α 2 s compared to result published in refs. [25,26]. We found that the second term proportional to s2F α 2 , whereas in refs. [25,26] it was listed as 32(1 − 4μ 2 )ζ ′ (−1)/ζ(−1). The author of refs. [25,26] has agreed that this was a typo in his article.

NNLO result -General case
It is relatively straightforward to generalize the previously obtained result (4.5) to the case that each quark has a separate chemical potential µ f . The final result is where the sums over f and g include all quark flavors, z f = 1/2 − iμ f , and Ω YM NNLO is the pure-glue contribution as before.

Mass prescription
As discussed in ref. [52], the two-loop perturbative electric gluon mass, first introduced by Braaten and Nieto in [22,23] is the most suitable for three-loop HTLpt calculations. We use the Braaten-Nieto (BN) mass prescription for m D in the remainder of the paper. Originally, the two-loop perturbative mass was calculated in refs. [22,23] for zero chemical potential, however, Vuorinen has generalized it to finite chemical potential. The resulting expression for m 2 D is [25,26] The effect of the in-medium quark mass parameter m q in thermodynamic functions is small and following ref. [52] we take m q = 0 which is the three loop variational solution.

Thermodynamic functions
In this section we present our final results for the NNLO HTLpt pressure, energy density, entropy density, trace anomaly, and speed of sound.

Running coupling
Below we will generally use the self-consistent one-loop running coupled implied by eq. (2.7), however, in some places we will try to gauge the sensitivity of the result to the order of the running coupling by comparing the impact of using one-or three-loop running. The three-loop running coupling can be expressed approximately as [130,131]  with t = ln(Λ 2 /Λ 2 MS ) and For one-loop running, we take b 1 = b 2 = 0. For both one-and three-loop running we fix the scale Λ MS by requiring that α s (1.5 GeV) = 0.326 which is obtained from lattice measurements [132]. For one-loop running, this procedure gives Λ MS = 176 MeV, and for three-loop running, one obtains Λ MS = 316 MeV.

Scales
For the renormalization scale we use separate scales, Λ g and Λ q , for purely-gluonic and fermionic graphs, respectively. We take the central values of these renormalization scales to be Λ g = 2πT and Λ = Λ q = 2π T 2 + µ 2 /π 2 . In all plots the thick lines indicate the result obtained using these central values and the light-blue band indicates the variation of the result under variation of both of these scales by a factor of two, e.g. πT ≤ Λ g ≤ 4πT .
For all numerical results below we use c A = N c = 3 and N f = 3.

Pressure
The QGP pressure can be obtained directly from the thermodynamic potential (4.5) where Λ above is understood to include both scales Λ g and Λ q . In  with lattice data from refs. [1,3,12]. In order to gauge the sensitivity of the results to the order of the running coupling, in fig. 1 we show the results obtained using a one-loop running and in fig. 2 the results obtained using a three-loop running. As can be seen by comparing these two sets, the sensitivity of the results to the order of the running coupling is small for T 250 MeV. As a result, unless the order of the running coupling turns out to have a significant effect on a given observable (see e.g. the fourth-order baryon number susceptibility), we will show the results obtained using a one-loop running coupling consistent with the counterterms necessary to renormalize the NNLO thermodynamic potential (2.7).
For an additional comparison we can compute the change in the pressure ∆P = P(T, Λ, µ) − P(T, Λ, 0) .  In figure 3 we plot ∆P as a function of the temperature for µ B = 300 MeV and µ B = 400 MeV. The solid lines are the NNLO HTLpt result and the dashed lines are the result obtained in the Stefan-Boltzmann limit. We note that in fig. 3 the lattice data from the Wuppertal-Budapest group [3] is computed up to O(µ 2 B ), whereas the HTLpt result includes all orders in µ B . As can be seen from this figure, the NNLO HTLpt result is quite close to the result obtained in the Stefan-Boltzmann limit. Note that the small correction in going from the Stefan-Boltzmann limit to NNLO HTLpt indicates that the fermionic sector is, to good approximation, weakly coupled for T 300 MeV.

Energy density
Once the pressure is known, it is straightforward to compute other thermodynamic functions such as the energy density by computing derivatives of the pressure with respect to the temperature and chemical potential. The energy density can be obtained via In figure 4 we plot the scaled NNLO HTLpt energy density for µ B = 0 (left) and µ B = 400 MeV (right) together with lattice data from refs. [1] and [4]. As we can see from this figure, there is reasonable agreement between the NNLO HTLpt energy density and the lattice data when the central value of the scale is used.

Entropy density
Similarly, we can compute the entropy density  We note that in the ideal gas limit, the entropy density becomes MeV lattice data are from [4]. For the HTLpt results a one-loop running coupling constant was used.

Trace anomaly
Since it is typically the trace anomaly itself which is computed on the lattice and then integrated to obtain the other thermodynamic functions, it is interesting to compare directly with lattice data for the trace anomaly. The trace anomaly is simply I = E − 3P. In the ideal gas limit, the trace anomaly goes to zero since E = 3P . When interactions are included, however, the trace anomaly (interaction measure) becomes non-zero. In figure  6 we plot the scaled NNLO HTLpt trace anomaly for µ B = 0 (left) and µ B = 400 MeV (right) together with lattice data from refs. [1] and [4]. As we can see from this figure, there is quite good agreement between the NNLO HTLpt trace anomaly and the lattice data for T 220 MeV when the central value of the scale is used.

Speed of sound
Another quantity which is phenomenologically interesting is the speed of sound. The speed of sound is defined as

Quark number susceptibilities
Having the full thermodynamic potential as a function of chemical potential(s) and temperature allows us to compute the quark number susceptibilities. In general, one can  Figure 8. The scaled second order baryon number susceptibility compared with various lattice data using one loop running (left) and three-loop running (right). The lattice data labeled WB, BNL-BI(B), BNL-BI(u,s), MILC, and TIFR come from refs. [2], [8], [9], [11], and [133], respectively. introduce a separate chemical potential for each quark flavor giving a N f -dimensional vector µ ≡ (µ 1 , µ 2 , ..., µ N f ). By taking derivatives of the pressure with respect to chemical potentials in this set, we obtain the quark number susceptibilities 5 Below we will use a shorthand notation for the susceptibilities by specifying derivatives by a string of quark flavors in superscript form, e.g. χ uu 2 = χ 200 , χ ds 2 = χ 011 , χ uudd 4 = χ 220 , etc.
When computing the derivatives with respect to the chemical potentials we treat Λ q as being a constant and only put the chemical potential dependence of the Λ q in after the derivatives are taken. We have done this in order to more closely match the procedure used to compute the susceptibilities using resummed dimensional reduction [43]. 6

Baryon number susceptibilities
We begin by considering the baryon number susceptibilities. The n th -order baryon number susceptibility is defined as  Figure 9. The scaled fourth order baryon number susceptibility compared with various lattice data using one loop running (left) and three-loop running (right). The lattice data labeled WB, BNL-BI(B), BNL-BI(u,s), MILC, and TIFR come from refs. [2], [8], [9], [11], and [133], respectively. For a three flavor system consisting of (u, d, s), the baryon number susceptibilities can be related to the quark number susceptibilities [14]  If we treat all quarks as having the same chemical potential µ u = µ d = µ s = µ = 1 3 µ B , eqs.  Figure 11. The N f = 2 + 1 NNLO HTLpt scaled sixth-order baryon susceptibility as a function of temperature.
wardly compute the baryon number susceptibility by computing derivatives of (4.5) with respect to µ.
In figure 8 we compare the NNLO HTLpt result for the second order baryon number susceptibility with lattice data from various groups. In the left panel of this figure we used the one-loop running and on the right we used the three-loop running. As one can see, for this quantity, the size of the light-blue band becomes larger if one uses the three-loop running, however, the central value obtained is very close in both cases.
Comparing to the lattice data we see that the NNLO HTLpt prediction is approximately 10% higher than the lattice data at T = 250 MeV and approximately 2% higher at T = 800 MeV. We note in this context that recently the four-loop second-order baryon number susceptibility has been computed in ref. [43] using the resummed dimensional reduction method. The result from this approach lies within the NNLO HTLpt scale variation band and is even closer to the lattice data with the error at T = 250 MeV being approximately 2% and 1% at T = 800 MeV. Our result, taken together with the resummed dimensional reduction results seem to indicate that the quark sector of the QGP can be quite accurately described using resummed perturbation theory for temperatures above approximately 300 MeV.
In figure 9 we compare the NNLO HTLpt result for the fourth order baryon number susceptibility with lattice data. Once again we show in the left and right panels, the result obtain using the one-loop running coupling and three-loop running coupling, respectively. Both the one-and three-loop running results are consistent with the lattice data shown; however, the lattice error bars on this quantity are somewhat large and the data are restricted to temperatures below 400 MeV, making it difficult to draw firm conclusions from this comparison. That being said, HTLpt makes a clear prediction for the temperature dependence of the fourth order baryon number susceptibility. It will be very interesting to see if future lattice data agree with this prediction.
susceptibilities as a function of temperature along with lattice data for this ratio. As we can see from this figure, this ratio very rapidly approaches the Stefan-Boltzmann limit if one considers the central NNLO HTLpt line. Comparing with the lattice data we see that the NNLO HTLpt result is below the lattice data for temperatures less than approximately 300 MeV. Without lattice data at higher temperatures, it's hard to draw a firm conclusion regarding the temperature at which HTLpt provides a good description of this quantity.
In figure 11 we show the NNLO HTLpt prediction for the sixth order baryon number susceptibility. To the best of our knowledge there is currently no publicly available lattice data for this quantity. It will be very interesting to see if these NNLO HTLpt predictions agree with lattice data as they becomes available.

Single quark number susceptibilities
We now consider the single quark number susceptibilities (7.1). For these we use the general expression for the NNLO thermodynamic potential with different quark chemical potentials (4.7). The resulting susceptibilities can either be diagonal (same flavor on all derivatives) or off-diagonal (different flavor on some or all indices). In HTLpt there are off-diagonal susceptibilities emerging explicitly from graphs F f 3c and F f 3j ; however, the latter vanishes when we use the variational mass prescription for the quark mass (m q = 0), so we need only consider the F f 3c graph. Additionally, there are potential off-diagonal contributions coming from all HTL terms since the Debye mass receives contributions from all quark flavors. In practice, however, because we evaluate derivatives with respect to the various chemical potentials and then take µ i → 0, one finds that all off-diagonal second order susceptibilities  Figure 13. Comparison of the N f = 2 + 1 NNLO HTLpt ratio of the fourth to second order single quark susceptibility with lattice data. For the HTLpt results a one-loop running coupling constant was used. The data labeled WB come from refs. [5,6]. vanish in HTLpt. Therefore, for the three-flavor case one has and, as a result, the single quark second order susceptibility is proportional to the baryon number susceptibility As a consequence, one can compute χ uuuu 4 directly from (4.7) or by computing χ B 4 using (4.5) and χ uudd 4 using (4.7) and applying the above relation. In our final plots we compute χ uuuu 4 directly from (4.7), however, we have checked that we obtain the same result if we use (7.7) instead.
In figure 12 (left) we plot our result for the fourth order single quark susceptibility χ uuuu 4 compared to lattice data from refs. [9], [8], [10], and [133]. As we can see from this figure, for the fourth order susceptibility there is very good agreement with available lattice data. In addition, the scale variation of the HTLpt result is quite small for this particular quantity. In figure 12 (right) we plot our result for the fourth order off-diagonal single quark susceptibility χ uudd 4 compared to lattice data. From this right panel we also see reasonably good agreement between the NNLO HTLpt result and the available lattice data.  Figure 14. The N f = 2 + 1 NNLO HTLpt scaled sixth-order diagonal and off-diagonal single quark susceptibilities χ 600 (left), χ 420 (middle), and χ 222 (right) as a function of temperature.
In figure 13 we plot the scaled ratio of the fourth-and second-order single quark susceptibilities. Once again we see good agreement between the NNLO HTLpt result and lattice data. Once again, for both figures 12 and 13, the lattice data are confined to relatively low temperatures. It will be interesting to compare higher temperature lattice data with the NNLO HTLpt prediction as they become available.
Finally, in figure 14 we plot the diagonal and off-diagonal sixth-order quark number susceptibilities χ 600 (left), χ 420 (middle), and χ 222 (right). We are unaware of any lattice data for these quantities. As before, it will be very interesting to see if these NNLO HTLpt predictions agree with lattice data as they becomes available.

Conclusions and outlook
In this paper, we presented the results of a NNLO (three-loop) HTLpt calculation of the thermodynamic potential of QCD at finite temperature and chemical potential(s). Our final result (4.7) is completely analytic and should be valid in the region of the phase diagram for which µ i 2πT . Based on the resulting thermodynamic potential we proceeded to calculate the pressure, energy density, entropy density, trace anomaly, and speed of sound of the QGP. In all cases we found very good agreement between the results obtained using the central values of the renormalization scales and available lattice data. Additionally, we have made predictions for the diagonal and off-diagonal sixth-order baryon number and single quark susceptibilities.
Looking to the future there are still many avenues for improvement in the HTLpt approach: (1) inclusion of the effects of finite quark masses (2) extension of results to µ i 2πT and eventually to T = 0, and (3) to potentially resum logarithms in order to reduce the scale variation of the final results (light-blue bands in all figures). Of these three, the second task is the most straightforward; however, in order to make more definitive and constrained statements it now seems necessary to start moving in directions (1) and (3) as well. In closing, we emphasize that HTLpt provides a gauge invariant reorganization of perturbation theory for calculating static quantities in thermal field theory. Since the NNLO HTLpt results are in good agreement with lattice data for various thermodynamic quantities down to temperatures that are relevant for LHC, it would therefore be interesting and challenging to apply HTLpt to the calculation of dynamic quantities, especially transport coefficients, at these temperatures.

A Expansion in mass parameters
In refs. [44,45] the NLO HTLpt thermodynamic potential was reduced to scalar sumintegrals. Evaluating these scalar sum-integrals exactly seems intractable, however, the sum-integrals can be calculated approximately by expanding them in powers of m D /T and m q /T following the method developed in ref. [32]. We will adopt the same strategy in this paper and include all terms through order g 5 assuming that m D and m q are O(g) at leading order. At each loop order, the contributions can be divided into those coming from hard and soft momenta, which are the momenta proportional to the scales T and gT respectively. In the one-loop diagrams, the contributions are either hard (h) or soft (s), while at the two-loop level, there are hard-hard (hh), hard-soft (hs), and soft-soft (ss) contributions. At three loops there are hard-hard-hard (hhh), hard-hard-soft (hhs), hard-soft-soft (hss), and soft-soft-soft (sss) contributions.

A.1 One-loop sum-integrals
We now review the mass expansion of the necessary one-loop sum-integrals considering separately the contributions from hard and soft momenta. We list the purely gluonic contributions when they are necessary for simpler exposition of the final result. Note that in order to simplify the results, when possible, it is best to add the corresponding iterated polarization and self-energy insertions that appear at higher order in δ, e.g. below we will also include F g 2d , F q 2d , F g 3m , and F f 3i as "one-loop" contributions.

Hard contributions
For one loop gluon (F g 1a ) and one loop ghost (F g 1b ) diagrams, we need to expand in order m 2 D : The one-loop graph with a gluon self-energy insertion (F g 2d ) has an explicit factor of m 2 D and, therefore, we only need to expand the sum-integral to first order in m 2 D : The one-loop graph with two gluon self-energy insertions (F g 3m ) must be expanded to zeroth order in m 2 The sum of eqs.
The one-loop fermionic graph F f 1b needs to expanded to second order in m 2 The one-loop fermion loop with a fermion self-energy insertion F f 2d must be expanded to first order in m 2 q , The one-loop fermion loop with two self-energy insertions F f 3i must be expanded to zeroth order in m 2 q : The sum of eqs. (A.5)-(A.7) is particularly simple This is the free energy of an ideal gas consisting of a single massless fermion.

Soft contributions
The soft contributions in the diagrams F 1a+1b , F 2d , and F 3m arise from the P 0 = 0 term in the sum-integral. 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 The total soft contribution from eqs. (A.9)-(A.11) is There are no soft contributions from the leading-order fermion diagrams or HTL counterterms (polarization and self-energy insertions).

A.2 Two-loop sum-integrals
For hard momenta, the self-energies are suppressed by m D /T and m q /T relative to the inverse free propagators, so we can expand in powers of Π T , Π L , and Σ. As was the case for the one-loop contributions, we once again treat the polarization and self-energy insertion NNLO diagrams as two-loop graphs in order to simplify the resulting expressions.

(hh) contribution
We first consider the contribution from fermionic diagrams. The (hh) contribution from F 2a and F 2b reads We consider next the (hh) contributions from F f 3d and F f 3f . The easiest way to calculate these term, is to expand the two-loop diagrams F f 2a and F f 2b to first order in m 2 D . This yields (A.14) Next we consider the (hh) contribution from the diagrams The sum of eqs.(A.13)-(A. 15) is For completeness, the hard-hard contribution coming from two-loop pure-glue diagrams is [49] F g(hh) α s π T 4 (A.17)

(hs) contribution
In the (hs) region, one gluon momentum is soft but the fermionic momentum is always hard. The terms that contribute through order g 2 m 3 D T and g 2 m 2 q m D T from F f 2a and F f 2b were calculated in ref. [44][45][46] and read The (hs) contribution from diagrams F f 3d and F f 3f can again be calculated from the diagrams F 2a and F 2b by Taylor expanding their contribution to first order in m 2 D . This yields We also need the (hs) contributions from the diagrams F f 3e , F f 3g , F f 3k and F f 3l . Again we calculate these contributions by expanding the two-loop diagrams F 2a and F 2b to first order in m 2 q . This yields The sum of eqs. .

A.3 Three-loop sum-integrals
We now list the mass-expanded sum-integrals necessary at three-loops. As before we organize the contributions according to whether the momentum flowing in a given propagator is hard or soft.

hhh contribution
The (hhh) contributions from diagrams F f 3a and F f 3b are Expressions for F f 3a , F f 3b , and F f 3c can be found in ref. [25,26]. The results are where ℵ(n, z) is defined in appendix D. The (hhh) contribution proportional to c A s F is where the integrals appearing above are evaluated in app. B.3. Finally, we note that there is no (hhh) contribution from F f 3o since this is a purely HTL diagram.

hhs contribution
The (hhs) contribution to the F 3a , F 3b , and F 3c+3j are Computing the necessary sum-integrals one finds Similarly, one obtains

hss contribution
The only three loop diagram involving a fermionic line that has a (hss) contribution is F f 3n which can be written as

B Sum-Integrals
We can define a set of "master" sum integral as in [25,26] (B.6)

B.1 One loop sum-integrals
The specific bosonic sun-integrals needed are The specific fermionic sun-integrals needed are and Using the two basic one-loop sum-integrals above, we can construct other one-loop sum-integrals that will be necessary here as follows: (B.12) This allows us to compute the following sum-integrals

B.2 Two loop sum-integrals
For the purposes of this paper we only need one new two-loop sum-integral

B.3 Three loop sum-integrals
The three-loop sum-integrals necessary are

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)

C.1 One-loop integrals
The general one-loop integral is given by (C.5)

C.2 Two-loop integrals
We also need a few two-loop integrals on the form J n = pq 1 p 2 + m 2 1 (q 2 + m 2 ) n 1 (p + q) 2 .
(C.6) Specifically, we need J 1 and J 2 which were calculated in refs. [49]: where n is assumed to be a non-negative integer and z is a general complex number given here by z = 1/2 − iμ. Above ζ denotes the Riemann zeta function, and Ψ is the digamma function Below we list Taylor expansions of the function ℵ and ℵ(n, z) for values of n necessary for calculation of the susceptibilities presented in the main text. For general application we evaluate the ℵ functions exactly using Mathematica.