text:
# Two-loop HTL Thermodynamics with Quarks

## I Introduction

The current generation of relativistic heavy-ion collision experiments should exceed the energy density necessary for the formation of a quark-gluon plasma. It is therefore necessary to have a quantitative theoretical framework which can be used to calculate the properties of a quark-gluon plasma. The usual line of reasoning is that since QCD is asymptotically free, its running coupling constant $`\alpha _s`$ becomes weaker as the temperature increases and therefore the behavior of hadronic matter at sufficiently high temperature should be calculable using perturbative methods. Unfortunately, a straightforward perturbative expansion in powers of $`\alpha _s`$ does not seem to be of any quantitative use even at temperatures many orders of magnitude higher than those achievable in heavy-ion collisions.

The problem can be seen by looking at the perturbative expansion of the free energy $``$ of a quark-gluon plasma, whose weak-coupling expansion has been calculated completely through order $`\alpha _s^{5/2}`$ AZ-95 ; KZ-96 ; BN-96

$``$ $`=`$ $`{\displaystyle \frac{8\pi ^2}{45}}T^4[_0+_2{\displaystyle \frac{\alpha _s}{\pi }}+_3\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^{3/2}`$ (1)
$`+_4\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^2+_5\left({\displaystyle \frac{\alpha _s}{\pi }}\right)^{5/2}+O\left(\alpha _s^3\mathrm{log}\alpha _s\right)],`$
with

$`_0`$ $`=`$ $`1+\frac{21}{32}N_f,`$ (2)
$`_2`$ $`=`$ $`{\displaystyle \frac{15}{4}}\left(1+\frac{5}{12}N_f\right),`$ (3)
$`_3`$ $`=`$ $`30\left(1+\frac{1}{6}N_f\right)^{3/2},`$ (4)
$`_4`$ $`=`$ $`237.2+15.97N_f0.413N_f^2`$ (5)
$`+{\displaystyle \frac{135}{2}}\left(1+\frac{1}{6}N_f\right)\mathrm{log}\left[{\displaystyle \frac{\alpha _s}{\pi }}\left(1+\frac{1}{6}N_f\right)\right]`$
$`{\displaystyle \frac{165}{8}}\left(1+{\displaystyle \frac{5}{12}}N_f\right)\left(1{\displaystyle \frac{2}{33}}N_f\right)\mathrm{log}{\displaystyle \frac{\mu }{2\pi T}},`$
$`_5`$ $`=`$ $`\left(1+\frac{1}{6}N_f\right)^{1/2}\left[799.2+21.96N_f+1.926N_f^2\right]`$ (6)
$`+{\displaystyle \frac{495}{2}}\left(1+{\displaystyle \frac{1}{6}}N_f\right)\left(1{\displaystyle \frac{2}{33}}N_f\right)\mathrm{log}{\displaystyle \frac{\mu }{2\pi T}},`$
where $`\mu `$ is the renormalization scale, $`\alpha _s=\alpha _s\left(\mu \right)`$ is the running coupling constant in the $`\overline{\text{MS}}`$ scheme, and we have set $`N_c=3`$. The coefficient of $`\alpha _s^3\mathrm{log}\alpha _s`$ has recently been computed KLRS-02 ; however, since there are unknown perturbative and non-perturbative contributions at $`O\left(\alpha _s^3\right)`$ we do not include terms higher than $`O\left(\alpha _s^{5/2}\right)`$ in Eq. (1).

In Fig. 1, the free energy with $`N_f=2`$ is shown as a function of the temperature $`T/T_c`$, where $`T_c`$ is the critical temperature for the deconfinement transition. In the plot we have scaled the free energy by the free energy of an ideal gas of quarks and gluons which for arbitrary $`N_c`$ and $`N_f`$ is

$$_{\mathrm{ideal}}=\frac{\pi ^2}{45}T^4\left(N_c^21+\frac{7}{4}N_cN_f\right).$$
(7)
The weak-coupling expansions through orders $`\alpha _s`$, $`\alpha _s^{3/2}`$, $`\alpha _s^2`$, and $`\alpha _s^{5/2}`$ are shown as bands that correspond to varying the renormalization scale, $`\mu `$, by a factor of two around the central value $`\mu =2\pi T`$. As successive terms in the weak-coupling expansion are added, the predictions change wildly and the sensitivity to the renormalization scale grows. It is clear that a reorganization of the perturbation series is essential if perturbative calculations are to be of any quantitative use at temperatures accessible in heavy-ion collisions.

The free energy can also be calculated nonperturbatively using lattice gauge theory Karsch . The thermodynamic functions for pure-glue QCD have been calculated with high precision by Boyd et al. lattice-0 . There have also been calculations which include dynamical quarks klp ; lattice-Nf . In Fig. 1 we have included the latest lattice estimate of Karsch et al. klp  for the free energy for $`N_f=2`$ flavors of light quarks. The band indicates the estimated systematic error of their result which is reported as (15$`\pm `$5)%. Note that the quarks in the simulations do have non-zero masses and that extrapolation to zero quark mass would require significant computing time. Due to the difficulty associated with the inclusion of light/massless dynamical quarks on the lattice it is therefore desirable to have analytic methods which can be used to estimate the thermodynamic functions.

The only rigorous method available for reorganizing perturbation theory in thermal QCD is dimensional reduction to an effective 3-dimensional field theory KLRS ; Kajantie-97 . The coefficients of the terms in the effective lagrangian are calculated using perturbation theory, but calculations within the effective field theory are carried out nonperturbatively using lattice gauge theory. Dimensional reduction has the same limitations as ordinary lattice gauge theory: it can be applied only to static quantities and only at zero baryon number density. Unlike in ordinary lattice gauge theory, light dynamical quarks do not require any additional computer power, because they only enter through the perturbatively calculated coefficients in the effective lagrangian. This method has been applied to the Debye screening mass for QCD Kajantie-97  as well as the pressure KLRS .

There are some proposals for reorganizing perturbation theory in QCD that are essentially just mathematical manipulations of the weak coupling expansion. The methods include Padé approximates Pade , Borel resummation Parwani , and self-similar approximates Yukalov . These methods are used to construct more stable sequences of successive approximations that agree with the weak-coupling expansion when expanded in powers of $`\alpha _s`$. These methods can only be applied to quantities for which several orders in the weak-coupling expansion are known, so they are limited in practice to the thermodynamic functions.

One promising approach for reorganizing perturbation theory in thermal QCD is to use a variational framework. The free energy $``$ is expressed as the variational minimum of a thermodynamic potential $`\mathrm{\Omega }(T,\alpha _s;m^2)`$ that depends on one or more variational parameters that we denote collectively by $`m^2`$:

$$(T,\alpha _s)=\mathrm{\Omega }(T,\alpha _s;m^2)|_{\mathrm{\Omega }/m^2=0}.$$
(8)
A particularly compelling variational formulation is the $`\mathrm{\Phi }`$-derivable approximation, in which the complete propagator is used as an infinite set of variational parameters Phi . The $`\mathrm{\Phi }`$-derivable thermodynamic potential $`\mathrm{\Omega }`$ is the 2PI effective action, the sum of all diagrams that are 2-particle-irreducible with respect to the complete propagator CJT-74 . The $`n`$-loop $`\mathrm{\Phi }`$-derivable approximations, in which $`\mathrm{\Omega }`$ is the the sum of 2PI diagrams with up to $`n`$ loops, form a systematically improvable sequence of variational approximations. Until recently, $`\mathrm{\Phi }`$-derivable approximations have proved to be intractable for relativistic field theories except for simple cases in which the self-energy is momentum-independent. However there has been some recent progress in solving the 3-loop $`\mathrm{\Phi }`$-derivable approximation for scalar field theories. Braaten and Petitgirard have developed an analytic method for solving the 3-loop $`\mathrm{\Phi }`$-derivable approximation for the massless $`\varphi ^4`$ field theory BP-01 . Van Hees and Knoll have developed numerical methods for solving the 3-loop $`\mathrm{\Phi }`$-derivable approximation for the massive $`\varphi ^4`$ field theory vHK-01 . They also investigated renormalization issues associated with the $`\mathrm{\Phi }`$-derivable approximation. These issues have recently been studied in detail by Blaizot, Iancu, and Reinosa BIR-03 .

The application of the $`\mathrm{\Phi }`$-derivable approximation to QCD was first discussed by McLerran and Freedman FM-77 . One problem with this approach is that the thermodynamic potential $`\mathrm{\Omega }`$ is gauge dependent, and so are the resulting thermodynamic functions. The gauge dependence is the same order in $`\alpha _s`$ as the truncation error when evaluated off the stationary point and twice the order in $`\alpha _s`$ when evaluated at the stationary point AS-02 . However the most serious problem is that even the application of 2-loop $`\mathrm{\Phi }`$-derivable approximation to gauge theories has proved to be intractable.

The 2-loop $`\mathrm{\Phi }`$-derivable approximation for QCD has been used as the starting point for HTL resummations of the entropy by Blaizot, Iancu and Rebhan BIR-99  and of the pressure by Peshier Peshier-00 . The thermodynamic potential $`\mathrm{\Omega }_{2\mathrm{loop}}`$ is a functional of the complete gluon propagator $`D_{\mu \nu }\left(P\right)`$. However, in order to make the problem tractable the authors in Refs. BIR-99  and Peshier-00  were forced to make a variational ansatz for the exact gluon propagator which they took as the HTL gluon propagator in the infrared and free in the ultraviolet with an aribitrary momentum scale separating the two momentum regions. Using this ansatz they were able to calculate the QCD thermodynamic functions; however, a first-principles calculation of the corrections to their results for gauge theories would require the inclusion of exact vertices as well as exact propagators thus making the problem intractable.

The difficulties in calculating quantities using $`\mathrm{\Phi }`$-derivable approximations in gauge theories motivates the use of simpler variational approximations. One such strategy that involves a single variational parameter $`m`$ has been called optimized perturbation theory Stevenson-81 , variational perturbation theory varpert , or the linear $`\delta `$ expansion deltaexp . This strategy was applied to the thermodynamics of the massless $`\varphi ^4`$ field theory by Karsch, Patkos and Petreczky under the name screened perturbation theory KPP-97 . The method has also been applied to spontaneously broken field theories at finite temperature CK-98 . The calculations of the thermodynamics of the massless $`\varphi ^4`$ field theory using screened perturbation theory have been extended to 3 loops ABS-01 . The calculations can be greatly simplified by using a double expansion in powers of the coupling constant and $`m/T`$ AS-01 .

HTL perturbation theory (HTLpt) is an adaptation of this strategy to thermal QCD ABS-99 . The exactly solvable theory used as the starting point is one whose propagators are the HTL quark and gluon propagators. The variational mass parameters $`m_D`$ and $`m_q`$ are identified with the Debye screening mass and the induced quark mass. The one-loop free energy in HTLpt was calculated for QCD in Ref. ABS-99  and for QCD with massless quarks in Ref. ABS-00 . At this order, the parameters $`m_D`$ and $`m_q`$ could not be determined variationally, so their perturbative limits were used. The resulting thermodynamic functions had errors of order $`\alpha _s`$, but the terms of order $`\alpha _s^{3/2}`$ associated with Debye screening were correct. A two-loop calculation is required to reduce the errors to order $`\alpha _s^2`$.

In a previous paper we calculated the thermodynamic functions of pure-glue QCD to next-to-leading order in HTLpt ABPS-02 . In that paper we showed that it was possible to renormalize the resulting expressions for the thermodynamic potential at next-to-leading order using only vacuum and mass counterterms and we also showed that the corrections to the thermodynamic functions in going from leading-order to next-to-leading order were small down to temperatures on the order of $`10T_c`$. In this paper we calculate the thermodynamic functions of QCD to next-to-leading order in HTLpt including the contributions from quark and quark-gluon interaction diagrams.

We begin with a brief summary of HTL perturbation theory including quarks in Section II. In Section III, we give the expressions for the one-loop and two-loop diagrams for the thermodynamic potential. In Section IV, we reduce those diagrams to scalar sum-integrals. We are unable to compute those sum-integrals exactly, so in Section V we evaluate them by treating $`m_D`$ and $`m_q`$ as $`O\left(g\right)`$ quantities and expand them in $`m_D/T`$ and $`m_q/T`$ keeping all terms that contribute up to $`O\left(g^5\right)`$. The diagrams are combined in Section VII to obtain the final result for the two-loop thermodynamic potential up to $`O\left(g^5\right)`$. In Section VIII, we present our numerical results for the free energy of QCD at leading and next-to-leading order in HTLpt. In Section IX we evaluate the free energy in the large $`N_f`$ limit where exact numerical results have been obtained moore ; ipp .

There are several appendices that contain technical details of the calculations. In Appendix A, we give the Feynman rules for HTL perturbation theory in Minkowski space to facilitate the application of this formalism to signatures of the quark-gluon plasma. The most difficult aspect of these calculations was the evaluation of the sum-integrals obtained from the expansion in $`m_D/T`$ and $`m_q/T`$. We give the results for these sum-integrals in Appendix B. The evaluation of some difficult thermal integrals that were required to obtain the sum-integrals is described in Appendix C.

## II HTL perturbation theory

The lagrangian density that generates the perturbative expansion for QCD can be expressed in the form

$`_{\mathrm{QCD}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{Tr}\left(G_{\mu \nu }G^{\mu \nu }\right)+i\overline{\psi }\gamma ^\mu D_\mu \psi `$ (9)
$`+_{\mathrm{gf}}+_{\mathrm{ghost}}+\mathrm{\Delta }_{\mathrm{QCD}}.`$
The gauge potential is $`A_\mu =A_\mu ^at^a`$, with generators $`t^a`$ of the fundamental representation of SU($`N_c`$) normalized so that $`\mathrm{Tr}t^at^b=\delta ^{ab}/2`$. The field strength tensor is $`G_{\mu \nu }=_\mu A_\nu _\nu A_\mu ig[A_\mu ,A_\nu ]`$. In the quark term there is an implicit sum over the $`N_f`$ quark flavors and $`D_\mu =_\mu igA_\mu `$ is the covariant derivative for the fundamental representation. The ghost term $`_{\mathrm{ghost}}`$ depends on the choice of the gauge-fixing term $`_{\mathrm{gf}}`$. Two choices for the gauge-fixing term that depend on an arbitrary gauge parameter $`\xi `$ are the general covariant gauge and the general Coulomb gauge:

$`_{\mathrm{gf}}`$ $`=`$ $`{\displaystyle \frac{1}{\xi }}\mathrm{Tr}\left(\left(^\mu A_\mu \right)^2\right)\text{covariant},`$ (10)
$`=`$ $`{\displaystyle \frac{1}{\xi }}\mathrm{Tr}\left(\left(𝐀\right)^2\right)\text{Coulomb}.`$ (11)
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 $`\left(F^a\right)^{bc}=if^{abc}`$ and generators of the fundamental representation as $`T^a`$ we define the following group theory factors:

$`\left[F^cF^c\right]^{ab}`$ $`=`$ $`f^{adc}f^{bdc}=c_A\delta ^{ab},`$
$`\mathrm{Tr}F^aF^b`$ $`=`$ $`s_A\delta ^{ab},`$
$`\delta ^{aa}`$ $`=`$ $`d_A,`$
$`\left[T^aT^a\right]_{ij}`$ $`=`$ $`c_F\delta _{ij},`$
$`\mathrm{Tr}T^aT^b`$ $`=`$ $`s_F\delta ^{ab},`$
$`\delta ^{ii}`$ $`=`$ $`d_F=s_Fd_A/c_F.`$ (12)
With the standard normalization

$`c_A=s_A`$ $`=`$ $`N_c,`$
$`d_A`$ $`=`$ $`N_c^21,`$
$`c_F`$ $`=`$ $`\left(N_c^21\right)/\left(2N_c\right),`$
$`s_F`$ $`=`$ $`N_f/2,`$
$`d_F`$ $`=`$ $`N_cN_f.`$ (13)
The perturbative expansion in powers of $`g`$ generates ultraviolet divergences. The renormalizability of perturbative QCD guarantees that all divergences in physical quantities can be removed by renormalization of the coupling constant $`\alpha _s=g^2/4\pi `$. There is no need for wavefunction renormalization, because physical quantities are independent of the normalization of the field. There is also no need for renormalization of the gauge parameter, because physical quantities are independent of the gauge parameter.

Hard-thermal-loop perturbation theory (HTLpt) is a reorganization of the perturbation series for thermal QCD. The lagrangian density is written as

$`=\left(_{\mathrm{QCD}}+_{\mathrm{HTL}}\right)|_{g\sqrt{\delta }g}+\mathrm{\Delta }_{\mathrm{HTL}}.`$ (14)
The HTL improvement term is

$`_{\mathrm{HTL}}={\displaystyle \frac{1}{2}}\left(1\delta \right)m_D^2\mathrm{Tr}\left(G_{\mu \alpha }{\displaystyle \frac{y^\alpha y^\beta }{\left(yD\right)^2}}_yG_\beta ^\mu \right)`$
$`+\left(1\delta \right)im_q^2\overline{\psi }\gamma ^\mu {\displaystyle \frac{y^\mu }{yD}}_y\psi ,`$ (15)
where in the first term $`D_\mu `$ is the covariant derivative in the adjoint representation, in the second term $`D_\mu `$ is the covariant derivative in the fundamental representation, $`y^\mu =(1,\widehat{𝐲})`$ is a light-like four-vector, and $`\mathrm{}_y`$ represents the average over the directions of $`\widehat{𝐲}`$. The term (II) has the form of the effective lagrangian that would be induced by a rotationally invariant ensemble of colored sources with infinitely high momentum. The parameter $`m_D`$ can be identified with the Debye screening mass and the parameter $`m_q`$ can be identified as the induced finite temperature quark mass. HTLpt is defined by treating $`\delta `$ as a formal expansion parameter.

The HTL perturbation expansion generates ultraviolet divergences. In QCD perturbation theory, renormalizability constrains the ultraviolet divergences to have a form that can be cancelled by the counterterm lagrangian $`\mathrm{\Delta }_{\mathrm{QCD}}`$. We will demonstrate that renormalized perturbation theory can be implemented by including a counterterm lagrangian $`\mathrm{\Delta }_{\mathrm{HTL}}`$ among the interaction terms in (14). There is no proof that the HTL perturbation expansion is renormalizable, so the general structure of the ultraviolet divergences is not known; however, it was shown in our previous paper ABPS-02  that it was possible to renormalize the next-to-leading order HTLpt prediction for the free energy of pure-glue QCD using only a vacuum counterterm and Debye mass counterterm. Here we show that when quarks are included it is also possible to renormalize the resulting expressions using only vacuum, Debye mass, and quark mass counterterms.

The leading term in the delta expansion of the vacuum energy, $`_0`$, counterterm $`\mathrm{\Delta }_0`$ was deduced in Ref. ABS-99  by calculating the free energy to leading order in $`\delta `$. The $`_0`$ counterterm $`\mathrm{\Delta }_0`$ must therefore have the form

$$\mathrm{\Delta }_0=\left(\frac{d_A}{128\pi ^2ϵ}+O\left(\delta \alpha _s\right)\right)\left(1\delta \right)^2m_D^4.$$
(16)
To calculate the free energy to next-to-leading order in $`\delta `$, we need the counterterm $`\mathrm{\Delta }_0`$ to order $`\delta `$ and the counterterms $`\mathrm{\Delta }m_D^2`$ and $`\mathrm{\Delta }m_q^2`$ to order $`\delta `$. We will show that there is a nontrivial cancellation of the ultraviolet divergences if the mass counterterms have the form

$`\mathrm{\Delta }m_D^2`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{3\pi ϵ}}\left[{\displaystyle \frac{11}{4}}c_As_F\right]m_D^2,`$ (17)
$`\mathrm{\Delta }m_q^2`$ $`=`$ $`{\displaystyle \frac{\alpha _s}{3\pi ϵ}}\left[{\displaystyle \frac{9}{8}}{\displaystyle \frac{d_A}{c_A}}\right]m_q^2.`$ (18)
Physical observables are calculated in HTLpt by expanding them in powers of $`\delta `$, truncating at some specified order, and then setting $`\delta =1`$. This defines a reorganization of the perturbation series in which the effects of the $`m_D^2`$ and $`m_q^2`$ terms in (II) are included to all orders but then systematically subtracted out at higher orders in perturbation theory by the $`\delta m_D^2`$ and $`\delta m_q^2`$ terms in (II). If we set $`\delta =1`$, the lagrangian (14) reduces to the QCD lagrangian (9). If the expansion in $`\delta `$ could be calculated to all orders, all dependence on $`m_D`$ and $`m_q`$ should disappear when we set $`\delta =1`$. However, any truncation of the expansion in $`\delta `$ produces results that depend on $`m_D`$ and $`m_q`$. Some prescription is required to determine $`m_D`$ and $`m_q`$ as a function of $`T`$ and $`\alpha _s`$. We choose to treat both as variational parameters that should be determined by minimizing the free energy. If we denote the free energy truncated at some order in $`\delta `$ by $`\mathrm{\Omega }(T,\alpha _s,m_D,m_q,\delta )`$, our prescription is

$`{\displaystyle \frac{}{m_D}}\mathrm{\Omega }\left(T,\alpha _s,m_D,m_q,\delta =1\right)`$ $`=`$ $`0,`$
$`{\displaystyle \frac{}{m_q}}\mathrm{\Omega }\left(T,\alpha _s,m_D,m_q,\delta =1\right)`$ $`=`$ $`0.`$ (19)
Since $`\mathrm{\Omega }\left(T,\alpha _s,m_D,m_q,\delta =1\right)`$ is a function of the variational parameters $`m_D`$ and $`m_q`$, we will refer to it as the thermodynamic potential. We will refer to the variational equations (19) as the gap equations. The free energy $``$ is obtained by evaluating the thermodynamic potential at the solution to the gap equations (19). Other thermodynamic functions can then be obtained by taking appropriate derivatives of $``$ with respect to $`T`$.

## III Diagrams for the thermodynamic potential

The thermodynamic potential at leading order in HTL perturbation theory for an $`SU\left(N_c\right)`$ gauge theory with $`N_f`$ massless quarks is

$`\mathrm{\Omega }_{\mathrm{LO}}=d_A_g+d_F_q+\mathrm{\Delta }_0_0,`$ (20)
where $`_g`$ is the contribution from each of the color states of the gluon:

$`_g`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _P}\left\{\left(d1\right)\mathrm{log}\left[\mathrm{\Delta }_T\left(P\right)\right]+\mathrm{log}\mathrm{\Delta }_L\left(P\right)\right\}.`$ (21)
The transverse and longitudinal HTL propagators $`\mathrm{\Delta }_T\left(P\right)`$ and $`\mathrm{\Delta }_L\left(P\right)`$ are given in (A.120) and (A.121). The quark contribution is

$`_q={\displaystyle _{\left\{P\right\}}}\mathrm{log}det[P/\mathrm{\Sigma }\left(P\right)],`$ (22)
where $`\mathrm{\Sigma }\left(P\right)`$ is the HTL fermion self-energy. The leading order counterterm $`\mathrm{\Delta }_0_0`$ was determined in Ref. ABS-99

$`\mathrm{\Delta }_0_0={\displaystyle \frac{d_A}{128\pi ^2ϵ}}m_D^4.`$ (23)
The thermodynamic potential at next-to-leading order in HTL perturbation theory can be written as

$`\mathrm{\Omega }_{\mathrm{NLO}}`$ $`=`$ $`\mathrm{\Omega }_{\mathrm{LO}}+d_A\left[_{3g}+_{4g}+_{gh}+_{gct}\right]`$ (24)
$`+d_As_F\left[_{3qg}+_{4qg}\right]+d_f_{qct}+\mathrm{\Delta }_1_0`$
$`+\mathrm{\Delta }_1m_D^2{\displaystyle \frac{}{m_D^2}}\mathrm{\Omega }_{\mathrm{LO}}+\mathrm{\Delta }_1m_q^2{\displaystyle \frac{}{m_q^2}}\mathrm{\Omega }_{\mathrm{LO}},`$
where $`\mathrm{\Delta }_1_0`$ and $`\mathrm{\Delta }_1m_D^2`$ are the terms of order $`\delta `$ in the vacuum energy density and mass counterterms. The contributions from the two-loop diagrams with the three-gluon and four-gluon vertices are

$`_{3g}`$ $`=`$ $`{\displaystyle \frac{c_A}{12}}g^2{\displaystyle _{PQ}}\mathrm{\Gamma }^{\mu \lambda \rho }(P,Q,R)`$ (25)
$`\times \mathrm{\Gamma }^{\nu \sigma \tau }(P,Q,R)\mathrm{\Delta }^{\mu \nu }\left(P\right)\mathrm{\Delta }^{\lambda \sigma }\left(Q\right)\mathrm{\Delta }^{\rho \tau }\left(R\right),`$
$`_{4g}`$ $`=`$ $`{\displaystyle \frac{c_A}{8}}g^2{\displaystyle _{PQ}}\mathrm{\Gamma }^{\mu \nu ,\lambda \sigma }(P,P,Q,Q)`$ (26)
$`\times \mathrm{\Delta }^{\mu \nu }\left(P\right)\mathrm{\Delta }^{\lambda \sigma }\left(Q\right),`$
where $`R=QP`$.

The contribution from the ghost diagram depends on the choice of gauge. The expressions in the covariant and Coulomb gauges are

$`_{gh}`$ $`=`$ $`{\displaystyle \frac{c_A}{2}}g^2{\displaystyle _{PQ}}{\displaystyle \frac{1}{Q^2}}{\displaystyle \frac{1}{R^2}}Q^\mu R^\nu \mathrm{\Delta }^{\mu \nu }\left(P\right)\text{covariant},`$ (28)
$`=`$ $`{\displaystyle \frac{c_A}{2}}g^2{\displaystyle _{PQ}}{\displaystyle \frac{1}{q^2}}{\displaystyle \frac{1}{r^2}}\left(Q^\mu Qnn^\mu \right)\left(R^\nu Rnn^\nu \right)`$
$`\times \mathrm{\Delta }^{\mu \nu }\left(P\right)\text{Coulomb}.`$
The contribution from the HTL gluon counterterm diagram is

$`_{\mathrm{gct}}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle _P}\mathrm{\Pi }^{\mu \nu }\left(P\right)\mathrm{\Delta }^{\mu \nu }\left(P\right).`$ (29)
The contributions from the two-loop diagrams with the quark-gluon three and four vertices are given by

$`_{3qg}`$ $`=`$ $`{\displaystyle \frac{1}{2}}g^2{\displaystyle _{P\left\{Q\right\}}}\text{Tr}[\mathrm{\Gamma }^\mu (P,Q,R)S\left(Q\right)`$ (30)
$`\times \mathrm{\Gamma }^\nu (P,Q,R)S\left(R\right)\left]\mathrm{\Delta }^{\mu \nu }\right(P)`$
$`_{4qg}`$ $`=`$ $`{\displaystyle \frac{1}{2}}g^2{\displaystyle _{P\left\{Q\right\}}}\text{Tr}\left[\mathrm{\Gamma }^{\mu \nu }(P,P,Q,Q)S\left(Q\right)\right]`$ (31)
$`\times \mathrm{\Delta }^{\mu \nu }\left(P\right),`$
where Tr implies taking the trace over $`\gamma `$-matrices. The contribution from the HTL quark counterterm is

$`_{\mathrm{qct}}={\displaystyle _{\left\{P\right\}}}\text{Tr}\left[\mathrm{\Sigma }\left(P\right)S\left(P\right)\right].`$ (32)
Provided that HTL perturbation theory is renormalizable, the ultraviolet divergences at any order in $`\delta `$ can be cancelled by renormalizations of the vacuum energy density $`_0`$, the HTL mass parameters $`m_D^2`$ and $`m_q^2`$, and the coupling constant $`\alpha _s`$. Renormalization of the coupling constant does not enter until order $`\delta ^2`$. We will calculate the thermodynamic potential as a double expansion in powers of $`m_D/T`$ and $`m_q/T`$, ad $`g`$ including all terms through 5<sup>th</sup> order. The $`\delta \alpha _s`$ term in $`\mathrm{\Delta }_0`$ does not contribute until 6<sup>th</sup> order in this expansion, so the term of order $`\delta `$ in $`\mathrm{\Delta }_0`$ can be obtained simply by expanding Eq. (23) to first order in $`\delta `$:

$`\mathrm{\Delta }_1_0={\displaystyle \frac{d_A}{64\pi ^2ϵ}}m_D^4.`$ (33)
The remaining ultraviolet divergences must be removed by renormalization of the mass parameters $`m_D`$ and $`m_q`$. We will show below that, at order $`\delta `$, all remaining divergences can be removed by the quark and Debye mass counterterms. This provides nontrivial evidence for the renormalizability of HTL perturbation theory at this order in $`\delta `$.

The sum of the 3-gluon, 4-gluon, and ghost contributions is gauge invariant. By using the Ward identities, one can easily show that the sum of these three diagrams is independent of the gauge parameter $`\xi `$. With more effort, one can show the equivalence of the covariant gauge expression with $`\xi =0`$ and the Coulomb gauge expression with $`\xi =0`$. In a similar manner, it can be shown that the sum of (30) and (31) is independent of $`\xi `$ within the class of covariant and Coloumb gauges, as well as the equivalence of the two with $`\xi =0`$.

## IV Reduction to scalar sum-integrals

The first step in calculating the quark contribution to the free energy is to reduce the sum of the diagrams to scalar-sum-integrals. The leading-order quark contribution can be rewritten as

$`_q=2{\displaystyle _{\left\{P\right\}}}\mathrm{log}P^22{\displaystyle _{\left\{P\right\}}}\mathrm{log}\left[{\displaystyle \frac{A_S^2A_0^2}{P^2}}\right],`$ (34)
where

$`A_0\left(P\right)`$ $`=`$ $`iP_0{\displaystyle \frac{m_q^2}{iP_0}}𝒯_P,`$ (35)
$`A_S\left(P\right)`$ $`=`$ $`\left|𝐩\right|+{\displaystyle \frac{m_q^2}{\left|𝐩\right|}}\left[1𝒯_P\right].`$ (36)
The HTL quark counterterm can be rewritten as

$`_{qct}=4{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{P^2+m_q^2}{A_S^2A_0^2}}.`$ (37)
We proceed to simplify the sum of (30) and (31) in Landau gauge. Using the Ward identities (A.114) and (A.117) the sum of (30) and (31) becomes

$`_{3qg+4qg}`$ $`=`$ $`{\displaystyle \frac{1}{2}}g^2{\displaystyle _{P\left\{Q\right\}}}\left\{\mathrm{\Delta }_X\left(P\right)\text{Tr}\left[\mathrm{\Gamma }^{00}S\left(Q\right)\right]\mathrm{\Delta }_T\left(P\right)\text{Tr}\left[\mathrm{\Gamma }^\mu S\left(Q\right)\mathrm{\Gamma }^\mu S\left(R^{}\right)\right]+\mathrm{\Delta }_X\left(P\right)\text{Tr}\left[\mathrm{\Gamma }^0S\left(Q\right)\mathrm{\Gamma }^0S\left(R^{}\right)\right]\right\},`$ (38)
where $`S`$ is the quark propagator, $`\mathrm{\Delta }_T`$ is the transverse gluon propagator, $`\mathrm{\Delta }_X`$ is a combinaiton of the transverse and longitudinal gluon propagators defined in (A.96), and $`R^{}=QP`$.

Performing the traces of $`\gamma `$-matrices gives

$`_{3qg+4qg}`$ $`=`$ $`g^2{\displaystyle _{P\left\{Q\right\}}}{\displaystyle \frac{1}{A_S^2\left(Q\right)A_0^2\left(Q\right)}}\{2(d1)\mathrm{\Delta }_T\left(P\right){\displaystyle \frac{\widehat{𝐪}\widehat{𝐫}A_S\left(Q\right)A_S\left(R\right)A_0\left(Q\right)A_0\left(R\right)}{A_S^2\left(R\right)A_0^2\left(R\right)}}`$ (39)
$`2\mathrm{\Delta }_X\left(P\right){\displaystyle \frac{A_0\left(Q\right)A_0\left(R\right)+A_S\left(Q\right)A_S\left(R\right)\widehat{𝐪}\widehat{𝐫}}{A_S^2\left(R\right)A_0^2\left(R\right)}}4m_q^2\mathrm{\Delta }_X\left(P\right){\displaystyle \frac{A_0\left(Q\right)A_s\left(Q\right)\widehat{𝐪}\widehat{𝐲}}{\left(PY\right)^2\left(QY\right)^2}}{\displaystyle \frac{1}{\left(QY\right)}}_{\widehat{𝐲}}`$
$`+{\displaystyle \frac{8m_q^2\mathrm{\Delta }_T\left(P\right)}{A_S^2\left(R\right)A_0^2\left(R\right)}}{\displaystyle \frac{\left(A_0\left(Q\right)A_S\left(Q\right)\widehat{𝐪}\widehat{𝐲}\right)\left(A_0\left(R\right)A_S\left(R\right)\widehat{𝐫}\widehat{𝐲}\right)}{\left(QY\right)\left(RY\right)}}_{\widehat{𝐲}}`$
$`+{\displaystyle \frac{4m_q^2\mathrm{\Delta }_X\left(P\right)}{A_S^2\left(R\right)A_0^2\left(R\right)}}{\displaystyle \frac{2A_0\left(R\right)A_S\left(Q\right)\widehat{𝐪}\widehat{𝐲}A_0\left(Q\right)A_0\left(R\right)A_S\left(Q\right)A_S\left(R\right)\widehat{𝐪}\widehat{𝐫}}{\left(QY\right)\left(RY\right)}}_{\widehat{𝐲}}\}+O(g^2m_q^4),`$
where $`A_0`$ and $`A_S`$ are defined in (A.109) and (A.110), respectively.

## V High-temperature Expansion

The free energy has been reduced to scalar sum-integrals. If we tried to evaluate the 2-loop HTL free energy exactly, there are terms that could at best be reduced to 5-dimensional integrals which would have to be evaluated numerically. We will therefore evaluate the sum-integrals approximately by expanding them in powers of $`m_D/T`$ and $`m_q/T`$. We will carry out the expansion to high enough order to include all terms through order $`g^5`$ if $`m_D`$ and $`m_q`$ are taken to be of order $`g`$.

The free energy can be divided into contributions from hard and soft momenta. We proceed to calculate the hard-hard and hard-soft contributions. There is no soft-soft contribution since one of the momenta in the loop is always fermionic and therefore hard.

### V.1 One-loop sum-integrals

The one-loop sum-integrals include the leading quark contribution (22) and the HTL quark counterterm (32). The leading order free energy must be expanded to order $`m_q^4`$ to include all terms through order $`g^5`$ if $`m_q`$ is taken to be of order $`g`$.

#### V.1.1 Hard contributions

The hard contribution from the LO gluon term (21) was given in ABPS-02  and reads

$`_g^{\left(h\right)}`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{45}}T^4+{\displaystyle \frac{1}{24}}\left[1+\left(2+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right)ϵ\right]`$ (40)
$`\times \left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}m_D^2T^2{\displaystyle \frac{1}{128\pi ^2}}\left({\displaystyle \frac{1}{ϵ}}7+2\gamma +{\displaystyle \frac{2\pi ^2}{3}}\right)`$
$`\times \left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}m_D^4.`$
The hard contribution from the HTL counterterm (29) was given in ABPS-02  and reads

$`_{gct}^{\left(h\right)}`$ $`=`$ $`{\displaystyle \frac{1}{24}}m_D^2T^2+{\displaystyle \frac{1}{64\pi ^2}}\left({\displaystyle \frac{1}{ϵ}}7+2\gamma +{\displaystyle \frac{2\pi ^2}{3}}\right)`$ (41)
$`\times \left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}m_D^4.`$
The sum-integrals over $`P`$ involve two momentum scales $`m_q`$ and $`T`$. Since $`P_0=\left(2n+1\right)\pi T`$, the momentum is always hard. We can therefore expand in powers of $`m_q^2`$. To second order in $`m_q^2`$, we obtain

$`_q^{\left(h\right)}`$ $`=`$ $`2{\displaystyle _{\left\{P\right\}}}\mathrm{log}P^24m_q^2{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{1}{P^2}}`$ (42)
$`+2m_q^4{\displaystyle _{\left\{P\right\}}}\left[{\displaystyle \frac{2}{P^4}}{\displaystyle \frac{1}{p^2P^2}}+{\displaystyle \frac{2}{p^2P^2}}𝒯_P{\displaystyle \frac{1}{p^2P_0^2}}\left(𝒯_P\right)^2\right].`$
Note that the function $`𝒯_P`$ cancels from the $`m_q^2`$ term. The values of the sum-integrals are given in Appendix B. Inserting those expressions, the hard quark contributions to the free energy reduce to

$`_q^{\left(h\right)}`$ $`=`$ $`{\displaystyle \frac{7\pi ^2}{180}}T^4+{\displaystyle \frac{1}{6}}\left[1+\left(22\mathrm{log}2+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right)ϵ\right]`$ (43)
$`\times \left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}m_q^2T^2+{\displaystyle \frac{1}{12\pi ^2}}\left(\pi ^26\right)m_q^4.`$
Note that this contribution is finite and so the leading order counterterm $`\mathrm{\Delta }_0_0`$ is the same as in the pure-glue case. The HTL quark counterterm is given in (37). Expanding this term to second order in $`m_q^2`$ yields

$`_{\mathrm{qct}}^{\left(h\right)}`$ $`=`$ $`4m_q^2{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{1}{P^2}}4m_q^4{\displaystyle _{\left\{P\right\}}}[{\displaystyle \frac{2}{P^4}}{\displaystyle \frac{1}{p^2P^2}}`$ (44)
$`+{\displaystyle \frac{2}{p^2P^2}}𝒯_P{\displaystyle \frac{1}{p^2P_0^2}}\left(𝒯_P\right)^2].`$
The values of the sum-integrals are given in Appendix B. Inserting those expressions, the hard contributions to the HTL quark counterterm reduce to

$`_{\mathrm{qct}}^{\left(h\right)}={\displaystyle \frac{1}{6}}m_q^2T^2{\displaystyle \frac{1}{6\pi ^2}}\left(\pi ^26\right)m_q^4.`$ (45)
Note that the first term in Eq. (45) cancels the order-$`ϵ^0`$ term in coefficient of $`m_q^2`$ in Eq. (43)

#### V.1.2 Soft contributions

The soft contribution comes from the $`P_0=0`$ term in the sum-integral. At soft momentum $`P=(0,𝐩)`$, the HTL self-energy functions reduce to $`\mathrm{\Pi }_T\left(P\right)=0`$ and $`\mathrm{\Pi }_L\left(P\right)=m_D^2`$. The transverse term vanishes in dimensional regularization because there is no momentum scale in the integral over $`𝐩`$. Thus the soft contribution comes from the longitudinal term only.

The soft contribution to the leading order free energy (21) was given in Ref. ABPS-02  and reads

$`_g^{\left(s\right)}`$ $`=`$ $`{\displaystyle \frac{1}{12\pi }}\left[1+{\displaystyle \frac{8}{3}}ϵ\right]\left({\displaystyle \frac{\mu }{2m_D}}\right)^{2ϵ}m_D^3T.`$ (46)
The soft contribution to the HTL gluon counterterm (29) given in Ref. ABPS-02  and reads

$`_{gct}^{\left(s\right)}`$ $`=`$ $`{\displaystyle \frac{1}{8\pi }}m_D^3T.`$ (47)
There are no soft contribution from the leading-order quark term (34) or from the HTL quark counterterm (37).

### V.2 Two-loop sum-integrals

The sum of the two-loop sum-integrals is given in (38). Since these integrals have an explicit factor of $`g^2`$, we need only to expand the sum-integrals to order $`m_q^2m_D/T^3`$ and $`m_D^3/T^3`$ to include all terms through order $`g^5`$.

The sum-integrals involve two momentum scales $`m_q,m_D`$ and $`T`$. In order to expand them in powers of.., we separate them into contributions from hard loop momenta and soft loop momenta. This gives two separate regions which we will denote $`\left(hh\right)`$ and $`\left(hs\right)`$. In the $`\left(hh\right)`$ region, all three momenta $`P`$, $`Q`$, and $`R`$ are hard. In the $`\left(hs\right)`$ region, two of the three momenta are hard and the other soft.

#### V.2.1 Contributions from the (hh) region

For hard momenta, the self-energies are suppressed by $`m_D^2/T^2`$ or $`m_q^2/T^2`$ relative to the propagators, so we can expand in powers of $`\mathrm{\Pi }_T`$, $`\mathrm{\Pi }_L`$, and $`\mathrm{\Sigma }`$.

The $`\left(hh\right)`$ contribution from (25)–(28) was given in Ref. ABPS-02  and reads

$`_{3g+4g+gh}^{\left(hh\right)}`$ $`=`$ $`{\displaystyle \frac{\pi ^2}{12}}{\displaystyle \frac{c_A\alpha _s}{3\pi }}T^4{\displaystyle \frac{7}{96}}\left[{\displaystyle \frac{1}{ϵ}}+4.621\right]{\displaystyle \frac{c_A\alpha _s}{3\pi }}`$ (48)
$`\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}m_D^2T^2.`$
The $`\left(hh\right)`$ contribution from (30) and (31) can be written as

$`_{3qg+4qg}^{\left(hh\right)}`$ $`=`$ $`\left(d1\right)g^2\left[{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{1}{P^2Q^2}}{\displaystyle _{P\left\{Q\right\}}}{\displaystyle \frac{2}{P^2Q^2}}\right]+2m_D^2g^2{\displaystyle _{P\left\{Q\right\}}}\left[{\displaystyle \frac{1}{p^2P^2Q^2}}𝒯_P+{\displaystyle \frac{1}{\left(P^2\right)^2Q^2}}{\displaystyle \frac{d2}{d1}}{\displaystyle \frac{1}{p^2P^2Q^2}}\right]`$ (49)
$`+m_D^2g^2{\displaystyle _{\left\{PQ\right\}}}\left[{\displaystyle \frac{d+1}{d1}}{\displaystyle \frac{1}{P^2Q^2r^2}}{\displaystyle \frac{4d}{d1}}{\displaystyle \frac{q^2}{P^2Q^2r^4}}{\displaystyle \frac{2d}{d1}}{\displaystyle \frac{PQ}{P^2Q^2r^4}}\right]𝒯_R`$
$`+m_D^2g^2{\displaystyle _{\left\{PQ\right\}}}\left[{\displaystyle \frac{3d}{d1}}{\displaystyle \frac{1}{P^2Q^2R^2}}+{\displaystyle \frac{2d}{d1}}{\displaystyle \frac{PQ}{P^2Q^2r^4}}{\displaystyle \frac{d+2}{d1}}{\displaystyle \frac{1}{P^2Q^2r^2}}+{\displaystyle \frac{4d}{d1}}{\displaystyle \frac{q^2}{P^2Q^2r^4}}{\displaystyle \frac{4}{d1}}{\displaystyle \frac{q^2}{P^2Q^2r^2R^2}}\right]`$
$`+2m_q^2g^2\left(d1\right){\displaystyle _{\left\{PQ\right\}}}\left[{\displaystyle \frac{1}{P^2Q_0^2Q^2}}+{\displaystyle \frac{p^2r^2}{P^2q^2Q_0^2R^2}}\right]𝒯_Q+2m_q^2g^2\left(d1\right){\displaystyle _{P\left\{Q\right\}}}\left[{\displaystyle \frac{2}{P^2\left(Q^2\right)^2}}{\displaystyle \frac{1}{P^2Q_0^2Q^2}}𝒯_Q\right]`$
$`+2m_q^2g^2\left(d1\right){\displaystyle _{\left\{PQ\right\}}}\left[{\displaystyle \frac{d+3}{d1}}{\displaystyle \frac{1}{P^2Q^2R^2}}{\displaystyle \frac{2}{P^2\left(Q^2\right)^2}}+{\displaystyle \frac{r^2p^2}{q^2P^2Q^2R^2}}\right].`$
Using the expressions for the sum-integrals in Appendix B, this reduces to

$`_{3qg+4qg}`$ $`=`$ $`{\displaystyle \frac{5\pi ^2}{72}}{\displaystyle \frac{\alpha _s}{\pi }}T^4{\displaystyle \frac{1}{72}}\left[{\displaystyle \frac{1}{ϵ}}+1.2963\right]{\displaystyle \frac{\alpha _s}{\pi }}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}m_D^2T^2+{\displaystyle \frac{1}{8}}\left[{\displaystyle \frac{1}{ϵ}}+8.96751\right]{\displaystyle \frac{\alpha _s}{\pi }}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}m_q^2T^2.`$ (50)
#### V.2.2 The (hs) contribution

In the $`\left(hs\right)`$ region, the momentum $`P`$ is soft. The momenta $`Q`$ and $`R`$ are always hard. The function that multiplies the soft propagator $`\mathrm{\Delta }_T(0,𝐩)`$ or $`\mathrm{\Delta }_X(0,𝐩)`$ can be expanded in powers of the soft momentum $`𝐩`$. In the case of $`\mathrm{\Delta }_T(0,𝐩)`$, the resulting integrals over $`𝐩`$ have no scale and they vanish in dimensional regularization. The integration measure $`_𝐩`$ scales like $`m_D^3`$, the soft propagator $`\mathrm{\Delta }_X(0,𝐩)`$ scales like $`1/m_D^2`$, and every power of $`p`$ in the numerator scales like $`m_D`$.

The $`\left(hs\right)`$ contribution from (25)–(28) was given in Ref. ABPS-02  and reads

$`_{3g+4g+gh}^{\left(hs\right)}`$ $`=`$ $`{\displaystyle \frac{\pi }{2}}{\displaystyle \frac{c_A\alpha _s}{3\pi }}m_DT^3{\displaystyle \frac{11}{32\pi }}\left({\displaystyle \frac{1}{ϵ}}+{\displaystyle \frac{27}{11}}+2\gamma \right){\displaystyle \frac{c_A\alpha _s}{3\pi }}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}\left({\displaystyle \frac{\mu }{2m_D}}\right)^{2ϵ}m_D^3T.`$ (51)
The only terms that contribute through order $`g^2m_D^3T`$ and $`m_q^2m_Dg^3T`$ from (30) and (31) are

$`_{3qg+4qg}^{\left(hs\right)}`$ $`=`$ $`g^2T{\displaystyle _𝐩}{\displaystyle \frac{1}{p^2+m_D^2}}{\displaystyle _{\left\{Q\right\}}}\left[{\displaystyle \frac{2}{Q^2}}{\displaystyle \frac{4q^2}{\left(Q^2\right)^2}}\right]+2m_D^2g^2T{\displaystyle _𝐩}{\displaystyle \frac{1}{p^2+m_D^2}}{\displaystyle _{\left\{Q\right\}}}\left[{\displaystyle \frac{1}{\left(Q^2\right)^2}}\left(3+d\right){\displaystyle \frac{2}{d}}{\displaystyle \frac{q^2}{\left(Q^2\right)^3}}+{\displaystyle \frac{8}{d}}{\displaystyle \frac{q^4}{\left(Q^2\right)^4}}\right]`$ (52)
$`4m_q^2g^2T{\displaystyle _𝐩}{\displaystyle \frac{1}{p^2+m_D^2}}{\displaystyle _{\left\{Q\right\}}}\left[{\displaystyle \frac{3}{\left(Q^2\right)^2}}{\displaystyle \frac{4q^2}{\left(Q^2\right)^3}}{\displaystyle \frac{4}{\left(Q^2\right)^2}}𝒯_Q{\displaystyle \frac{2}{Q^2}}{\displaystyle \frac{1}{\left(QY\right)^2}}_{\widehat{𝐲}}\right].`$
In the terms that are already of order $`g^2m_D^3T`$, we can set $`R=Q`$. In the terms of order $`g^2m_DT^3`$, we must expand the sum-integrand to second order in $`𝐩`$. After averaging over angles of $`𝐩`$, the linear terms in $`𝐩`$ vanish and quadratic terms of the form $`p^ip^j`$ are replaced by $`p^2\delta ^{ij}/d`$. We can set $`p^2=m_D^2`$, because any factor proportional to $`p^2+m_D^2`$ will cancel the denominator of the integral over $`𝐩`$, leaving an integral with no scale. This gives

$`_{3qg+4qg}^{\left(hs\right)}`$ $`=`$ $`{\displaystyle \frac{1}{6}}\alpha _sm_DT^3+{\displaystyle \frac{\alpha _s}{24\pi ^2}}\left[{\displaystyle \frac{1}{ϵ}}+1+2\gamma +4\mathrm{log}2\right]`$
$`\times \left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}\left({\displaystyle \frac{\mu }{2m_D}}\right)^{2ϵ}m_D^3T{\displaystyle \frac{\alpha _s}{2\pi ^2}}m_q^2m_DT.`$
#### V.2.3 The $`(ss)`$ contributions

The $`\left(ss\right)`$ contribution from (25)–(28) was given in Ref. ABPS-02  and reads

$`_{3g+4g+gh}^{\left(ss\right)}`$ $`=`$ $`{\displaystyle \frac{3}{16}}\left[{\displaystyle \frac{1}{ϵ}}+3\right]{\displaystyle \frac{c_A\alpha _s}{3\pi }}\left({\displaystyle \frac{\mu }{2m_D}}\right)^{4ϵ}m_D^2T^2.`$ (54)
There is no $`\left(ss\right)`$ contribution from the diagrams involving fermions.

## VI HTL-Improved Thermodynamics

The free energy at second order in HTL perturbation theory defines a function $`\mathrm{\Omega }\left(T,\alpha _s,m_D,m_q,\delta =1\right)`$. We will refer to this function as the thermodynamic potential. To obtain the free energy $`F\left(T\right)`$ as a function of the temperature, we need to specify a prescription for the mass parameter $`m_D`$ as a function of $`T`$ and $`\alpha _s`$.

## VII Thermodynamic Potential

In this section, we calculate the thermodynamic potential $`\mathrm{\Omega }\left(T,\alpha _s,m_D,m_q,\delta =1\right)`$ explicitly, first to leading order in the $`\delta `$ expansion and then to next-to-leading order.

### VII.1 Leading order

The complete expression for the leading order thermodynamic potential is the sum of the contributions from 1-loop diagrams and the leading term (23) in the vacuum energy counterterm. The contributions from the 1-loop diagrams, including all terms through order $`g^5`$, is the sum of (40), (43), and (46)

$`\mathrm{\Omega }_{1\mathrm{loop}}`$ $`=`$ $`d_A{\displaystyle \frac{\pi ^2T^4}{45}}\{1+{\displaystyle \frac{7}{4}}{\displaystyle \frac{d_F}{d_A}}{\displaystyle \frac{15}{2}}\widehat{m}_D^230{\displaystyle \frac{d_F}{d_A}}\widehat{m}_q^2`$ (55)
$`+30\widehat{m}_D^3+{\displaystyle \frac{45}{8}}\left({\displaystyle \frac{1}{ϵ}}+2\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}7+2\gamma +{\displaystyle \frac{2\pi ^2}{3}}\right)\widehat{m}_D^4`$
$`60{\displaystyle \frac{d_F}{d_A}}(\pi ^26)\widehat{m}_q^4\},`$
where $`\widehat{m}_D`$, $`\widehat{m}_q`$ and $`\widehat{\mu }`$ are dimensionless variables:

$`\widehat{m}_D`$ $`=`$ $`{\displaystyle \frac{m_D}{2\pi T}},`$ (56)
$`\widehat{m}_q`$ $`=`$ $`{\displaystyle \frac{m_q}{2\pi T}},`$ (57)
$`\widehat{\mu }`$ $`=`$ $`{\displaystyle \frac{\mu }{2\pi T}}.`$ (58)
Adding the counterterm (23), we obtain the thermodynamic potential at leading order in the delta expansion:

$`\mathrm{\Omega }_{\mathrm{LO}}`$ $`=`$ $`d_A{\displaystyle \frac{\pi ^2T^4}{45}}\{1+{\displaystyle \frac{7}{4}}{\displaystyle \frac{d_F}{d_A}}{\displaystyle \frac{15}{2}}\widehat{m}_D^230{\displaystyle \frac{d_F}{d_A}}\widehat{m}_q^2`$ (59)
$`+30\widehat{m}_D^3+{\displaystyle \frac{45}{4}}\left(\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}{\displaystyle \frac{7}{2}}+\gamma +{\displaystyle \frac{\pi ^2}{3}}\right)\widehat{m}_D^4`$
$`60{\displaystyle \frac{d_F}{d_A}}(\pi ^26)\widehat{m}_q^4\},`$
### VII.2 Next-to-leading order

The complete expression for the next-to-leading order correction to the thermodynamic potential is the sum of the contributions from all 2-loop diagrams, the quark and gluon counterterms, and renormalization counterterms. The contributions from the 2-loop diagrams, including all terms though order $`g^5`$ is the sum of (48), (50), (51), (LABEL:F2loop-qhs), and (54) multiplied by the appropriate group structure constants listed in (24):

$`\mathrm{\Omega }_{2\mathrm{loop}}`$ $`=`$ $`d_A{\displaystyle \frac{\pi ^2T^4}{45}}{\displaystyle \frac{\alpha _s}{\pi }}\{{\displaystyle \frac{5}{4}}(c_A+{\displaystyle \frac{5}{2}}s_F)+15(c_A+s_F)\widehat{m}_D{\displaystyle \frac{55}{8}}[(c_A{\displaystyle \frac{4}{11}}s_F)({\displaystyle \frac{1}{ϵ}}+4\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}})`$ (60)
$`c_A({\displaystyle \frac{72}{11}}\mathrm{log}\widehat{m}_D1.96869)0.4714s_F]\widehat{m}_D^2{\displaystyle \frac{45}{2}}s_F[{\displaystyle \frac{1}{ϵ}}+4\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}+8.96751]\widehat{m}_q^2+180s_F\widehat{m}_D\widehat{m}_q^2`$
$`+{\displaystyle \frac{165}{4}}[(c_A{\displaystyle \frac{4}{11}}s_F)({\displaystyle \frac{1}{ϵ}}+4\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}2\mathrm{log}\widehat{m}_D)+c_A({\displaystyle \frac{27}{11}}+2\gamma ){\displaystyle \frac{4}{11}}s_F(1+2\gamma +4\mathrm{log}2)]\widehat{m}_D^3\}.`$
The HTL gluon counterterm is the sum of (41) and (47)

$`\mathrm{\Omega }_{\mathrm{gct}}`$ $`=`$ $`d_A{\displaystyle \frac{\pi ^2T^4}{45}}\left\{{\displaystyle \frac{15}{2}}\widehat{m}_D^245\widehat{m}_D^3{\displaystyle \frac{45}{4}}\left({\displaystyle \frac{1}{ϵ}}+2\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}7+2\gamma +{\displaystyle \frac{2\pi ^2}{3}}\right)\widehat{m}_D^4\right\}.`$ (61)
The HTL quark counterterm is given by (45)

$`\mathrm{\Omega }_{\mathrm{qct}}`$ $`=`$ $`d_F{\displaystyle \frac{\pi ^2T^4}{45}}\left\{30\widehat{m}_q^2+120\left(\pi ^26\right)\widehat{m}_q^4\right\}.`$ (62)
The ultraviolet divergences that remain after these 3 terms are added can be removed by renormalization of the vacuum energy density $`_0`$ and the HTL mass parameter $`m_D`$. The renormalization contributions at first order in $`\delta `$ are

$`\mathrm{\Delta }\mathrm{\Omega }=\mathrm{\Delta }_1_0+\mathrm{\Delta }_1m_D^2{\displaystyle \frac{}{m_D^2}}\mathrm{\Omega }_{\mathrm{LO}}+\mathrm{\Delta }_1m_q^2{\displaystyle \frac{}{m_q^2}}\mathrm{\Omega }_{\mathrm{LO}}.`$ (63)
Using the results listed in Eqs. (17), (18), and (33) the complete contribution from the counterterm at first order in $`\delta `$ is

$`\mathrm{\Delta }\mathrm{\Omega }`$ $`=`$ $`d_A{\displaystyle \frac{\pi ^2T^4}{45}}\{{\displaystyle \frac{45}{4ϵ}}\widehat{m}_D^4+{\displaystyle \frac{\alpha _s}{\pi }}[{\displaystyle \frac{55}{8}}(c_A{\displaystyle \frac{4}{11}}s_F)({\displaystyle \frac{1}{ϵ}}+2\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}+2)\widehat{m}_D^2`$
$`{\displaystyle \frac{165}{4}}(c_A{\displaystyle \frac{4}{11}}s_F)({\displaystyle \frac{1}{ϵ}}+2\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}2\mathrm{log}\widehat{m}_D+2)\widehat{m}_D^3+{\displaystyle \frac{45}{2}}s_F({\displaystyle \frac{1}{ϵ}}+2+2\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}2\mathrm{log}2+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}})\widehat{m}_q^2]\}.`$
Adding the contributions from the two-loop diagrams in (60), the HTL gluon and quark counterterms in (61) and (62), the contribution from vacuum and mass renormalizations in (LABEL:OmegaVMct), and the leading order thermodynamic potential in (59) we obtain the complete expression for the QCD thermodynamic potential at next-to-leading order in HTLpt:

$`\mathrm{\Omega }_{\mathrm{NLO}}`$ $`=`$ $`d_A{\displaystyle \frac{\pi ^2T^4}{45}}\{1+{\displaystyle \frac{7}{4}}{\displaystyle \frac{d_F}{d_A}}15\widehat{m}_D^3{\displaystyle \frac{45}{4}}(\mathrm{log}\widehat{{\displaystyle \frac{\mu }{2}}}{\displaystyle \frac{7}{2}}+\gamma +{\displaystyle \frac{\pi ^2}{3}})\widehat{m}_D^4+60{\displaystyle \frac{d_F}{d_A}}(\pi ^26)\widehat{m}_q^4`$ (65)
$`+{\displaystyle \frac{\alpha _s}{\pi }}[{\displaystyle \frac{5}{4}}(c_A+{\displaystyle \frac{5}{2}}s_F)+15(c_A+s_F)\widehat{m}_D{\displaystyle \frac{55}{4}}\{c_A(\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}{\displaystyle \frac{36}{11}}\mathrm{log}\widehat{m}_D2.001){\displaystyle \frac{4}{11}}s_F(\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}2.337)\}\widehat{m}_D^2`$
$`45s_F(\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}+2.192)\widehat{m}_q^2+{\displaystyle \frac{165}{2}}\{c_A(\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}+{\displaystyle \frac{5}{22}}+\gamma ){\displaystyle \frac{4}{11}}s_F(\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}{\displaystyle \frac{1}{2}}+\gamma +2\mathrm{log}2)\}\widehat{m}_D^3+180s_F\widehat{m}_D\widehat{m}_q^2]\}.`$
### VII.3 Gap Equation

The quark and gluon mass parameters, $`m_q`$ and $`m_D`$, are determined variationally by requiring that the derivative of $`\mathrm{\Omega }_{\mathrm{NLO}}`$ with respect to each parameter taken holding the other constant vanishes

$`{\displaystyle \frac{}{m_q}}\mathrm{\Omega }_{\mathrm{NLO}}(T,\alpha _s,m_D,m_q)`$ $`=`$ $`0,`$ (66)
$`{\displaystyle \frac{}{m_D}}\mathrm{\Omega }_{\mathrm{NLO}}(T,\alpha _s,m_D,m_q)`$ $`=`$ $`0.`$ (67)
The first equation above results in the following equation for $`m_q`$

$$8\frac{d_F}{d_A}\left(\pi ^26\right)\widehat{m}_q^2=\frac{\alpha _ss_F}{\pi }\left[3\left(\mathrm{log}\frac{\widehat{\mu }}{2}+2.192\right)12\widehat{m}_D\right].$$
(68)
In the limit of small $`\alpha _s`$ the above gap equation does not go to the perturbative limit for the quark mass which is $`\widehat{m}_{q,\mathrm{pert}}^2=c_F\alpha _s/8\pi `$. The fact that $`m_q`$ does not go to the perturbative value in the small $`\alpha _s`$ limit is due to the fact that the perturbative limit of the quark gap equation results from terms which are $`O\left(\alpha _s^2\right)`$ and these terms are not included completely at NLO in HTLpt. One might hope that going to the next order in HTLpt would cure this problem; however, this will in fact not happen since the fermion sector is infrared safe and therefore only even powers of $`\widehat{m}_q`$ will appear at each order. At NNLO all terms contributing at $`O\left(\alpha _s^2\right)`$ at NLO will be replaced by explicit powers of $`\alpha _s`$ and all $`\widehat{m}_q`$ dependence will be pushed up to $`O\left(\alpha _s^3\right)`$. This behavior will persist at all orders in HTLpt so that at any order the weak-coupling limit of the gap equation quark mass will be scale dependent. In order to circumvent this problem we can consider other possible prescriptions for choosing $`\widehat{m}_q`$ which include requiring that $`\widehat{m}_q`$ be equal to its perturbative value for all $`\alpha _s`$ or requiring that $`\widehat{m}_q`$ be proportional to $`\widehat{m}_D`$ with the proportionality constant fixed in the weak-coupling limit.

Performing the derivative with respect to $`m_D`$ while holding $`m_q`$ fixed results in the following gap equation for $`m_D`$

$`45\widehat{m}_D^2\left[1+\left(\mathrm{log}\widehat{{\displaystyle \frac{\mu }{2}}}{\displaystyle \frac{7}{2}}+\gamma +{\displaystyle \frac{\pi ^2}{3}}\right)\widehat{m}_D\right]`$
$`={\displaystyle \frac{\alpha }{\pi }}\{15(c_A+s_F){\displaystyle \frac{55}{2}}[c_A(\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}{\displaystyle \frac{36}{11}}\mathrm{log}\widehat{m}_D3.637)`$
$`{\displaystyle \frac{4}{11}}s_F(\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}2.337)]\widehat{m}_D+{\displaystyle \frac{495}{2}}[c_A(\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}+{\displaystyle \frac{5}{22}}+\gamma )`$
$`{\displaystyle \frac{4}{11}}s_F(\mathrm{log}{\displaystyle \frac{\widehat{\mu }}{2}}{\displaystyle \frac{1}{2}}+\gamma +2\mathrm{log}2)]\widehat{m}_D^2+180s_F\widehat{m}_q^2\}.`$ (69)
The last term in Eq. (69) proportional to $`\widehat{m}_q^2`$ can be written in terms of $`\widehat{m}_D`$ using Eq. (68). In Fig. 3 we plot the solutions to the gap equations for $`m_D`$ and $`m_q`$ for $`N_c=3`$ and $`N_f=2`$. The solution for $`m_D`$ goes to the perturbative value in the limit of small $`\alpha _s`$, decreases below the perturbative value as $`\alpha _s`$ increases, and becomes larger than the perturbative value at $`\alpha _s0.11`$. The solution for $`m_q`$ does not go to the perturbative value in the limit of small $`\alpha _s`$ and is instead scale dependent even at lowest order as dicussed above. As $`\alpha _s`$ increases the value of $`m_q`$ remains very flat regardless of the scale, changing significantly only near $`\alpha _s0.10`$.

## VIII Free energy

The free energy is obtained by evaluating the leading and next-to-leading order thermodynamic potentials, (59) and (65), at the solution to the gap equations (68) and (69). In Fig. 4 we plot the leading and next-to-leading order HTLpt predictions for the free energy of QCD with $`N_c=3`$ and $`N_f=2`$. We have studied the alternative prescriptions for the quark mass discussed in the previous section and find that the NLO free energy obtained using these prescriptions is numerically indistinguishable from that obtained using the quark gap equations. As can be seen from this figure the corrections in going from LO to NLO are small over the entire temperature range, especially when compared to convergence of the perturbative result. Additionally, the scale variation of the NLO HTLpt result for the free energy is much smaller than the LO showing that the partial resummation of the scale dependent logarithms reduces the scale variation of the final results significantly.

However, as was the case in pure-glue QCD ABPS-02 , the results seem to lie significantly above the lattice data which is available below $`5T_c`$. There are several reasons for why HTLpt might fail to describe the lattice data in this temperature range. One is that the hard modes are not resummed properly within HTLpt and that a description using a $`\mathrm{\Phi }`$-derivable approach which explicitly separates the hard and soft modes as done in Ref. BIR-99  is better. A second is that HTLpt discards some important physics like topological modes or the $`Z_N`$ symmetry of QCD near the phase transition.

A third possibility is that the expansion in $`m_D/T`$ and $`m_q/T`$ breaks down at these temperatures. Numerically, $`m_D/T1.2`$ and $`m_q/T0.5`$ at $`5T_c`$ and $`m_D/T1.6`$ and $`m_q/T0.6`$ at $`2T_c`$, which casts doubt on the applicability of the expansion in this temperature range. However, in the case of pure-glue we have been able to compare the LO HTLpt result expanded to $`O\left(\widehat{m}_D^6\right)`$ with the non-truncated LO expression which is accurate to all orders in $`\widehat{m}_D`$ and find that the expansions converge very rapidly. Numerically we find that at $`m_D/T=5`$ truncations of the LO order result accurate to $`\widehat{m}_D^4`$ and $`\widehat{m}_D^6`$ reproduce the exact result to 5% and 0.2%, respectively. There have also been studies of the convergence of the mass expansions of the three-loop free energy for a massless scalar field theory using screened perturbation theory AS-01  and the $`\mathrm{\Phi }`$-derivable approach BP-01  which demonstrated that mass expansions also converge very rapidly at NLO and NNLO. This gives us some confidence that the truncated NLO solutions are numerically reliable.

## IX Large $`N_f`$

In the limit that $`N_f`$ is taken large while holding $`g^2N_f`$ fixed it is possible to solve for the $`O\left(N_f^0\right)`$ contribution to the free energy exactly moore ; ipp . In Fig. 5 we plot the NLO HTLpt prediction for the $`O\left(N_f^0\right)`$ contribution to the free energy along with the numerical result of Ref. ipp  and the perturbative prediction which is accurate to $`O\left(\alpha _s^{5/2}\right)`$. In Fig. 6 we plot the NLO HTLpt prediction for $`m_D`$ in the large $`N_f`$ limit and the exact numerical result tonyprivate . The HTLpt predictions for both the free energy and the Debye mass seem diverge from the exact result around $`g_{\mathrm{eff}}2`$ regardless of the scale which is chosen; however, for both quantities, choosing the scale to be $`\mu =\mu _{\mathrm{DR}}=\pi e^\gamma T`$ seems to provide a reasonable reproduction of the exact results. This result is comparable to the performance of the $`\mathrm{\Phi }`$-derivable approach in the large $`N_f`$ limit rebhan-03 .

## X Conclusions

In this paper we have extended our previous HTLpt calculation of the thermodynamic functions in pure-glue QCD to include the contribution of $`N_f`$ massless quarks. We have presented results for the leading- and next-to-leading-order HTLpt predictions for the QCD free energy for arbitrary $`N_f`$. Using the NLO HTLpt expression for the thermodynamic potential we were able to find variational solutions for both the quark and gluon mass parameters allowing a first-principles prediction of the QCD free energy. As in the case of pure-glue we find that the NLO HTLpt prediction lies significantly above the available lattice data below $`5T_c`$; however, the problem of oscillation of successive approximations and large scale dependence of the perturbative results is eliminated by using this reorganization.

The failure of HTLpt to describe the lattice data in this region could be attributed to the failure of the expansion performed in $`\widehat{m}_D`$ and $`\widehat{m}_q`$; however, a study of the convergence of the truncated LO expressions to a numerical evaluation of the exact LO expression shows that these expansions converge very rapidly. Therefore, we are steered towards the conclusion that a systematic description of QCD thermodynamics using HTLpt is not appropriate below $`5T_c`$. The $`\mathrm{\Phi }`$-derivable approach seems to agree better with the lattice data in this range so perhaps HTLpt is not resumming the hard modes properly and an explicit separation of hard and soft scales is required. However, we should point out that some authors believe that a description of QCD thermodynamics near the phase transition in terms of Polyakov loops is necessary pisarski-00 .

We have also compared the NLO HTLpt free energy and Debye mass with exact results which are available in the large $`N_f`$ limit. This comparison shows that, in the large $`N_f`$ limit, NLO HTLpt agrees with the exact result only out to $`g_{\mathrm{eff}}2`$ and has large scale dependence after this point. The large scale dependence is not surprising given the fact that in the large $`N_f`$ limit the running of the coupling constant is driven by the presence of the Landau singularity and even the exact results are sensitive to this beyond $`g_{\mathrm{eff}}5`$. The poor performance of NLO HTLpt, however, is comparable to recent large $`N_f`$ predictions within the $`\mathrm{\Phi }`$-derivable approach. The failure of both approaches to agree better with the exact result for large values of $`g_{\mathrm{eff}}`$ is an indication that a description of strongly-coupled QCD thermodynamics solely in terms of HTL quasiparticles is perhaps inappropriate. However, it is possible that the physics of large-$`N_f`$ QCD is so different from that of QCD with a small number of flavors that it cannot serve as a definitive testing ground for the applicability of the quasiparticle approach to the physical case peshier-03 .

## Acknowledgments

M.S. would like to thank A.Rebhan for discussions and comments. J.O.A was supported by the Stichting voor Fundamenteel Onderzoek der Materie (FOM), which is supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO). M.S. was supported by US DOE grant DE-FG02-96ER40945 and by an FWF Der Wissenschaftsfonds Lise Meitner fellowship M689.

## Appendix A HTL Feynman Rules

In this appendix, we present Feynman rules for HTL perturbation theory in QCD. We give explicit expressions for the propagators and for the quark-gluon 3 and 4 vertices. The Feynman rules are given in Minkowski space to facilitate applications to real-time processes. A Minkowski momentum is denoted $`p=(p_0,𝐩)`$, and the inner product is $`pq=p_0q_0𝐩𝐪`$. The vector that specifies the thermal rest frame is $`n=(1,\mathrm{𝟎})`$.

### A.1 Gluon Self-energy

The HTL gluon self-energy tensor for a gluon of momentum $`p`$ is

$`\mathrm{\Pi }^{\mu \nu }\left(p\right)=m_D^2\left[𝒯^{\mu \nu }(p,p)n^\mu n^\nu \right].`$ (A.70)
The tensor $`𝒯^{\mu \nu }(p,q)`$, which is defined only for momenta that satisfy $`p+q=0`$, is

$`𝒯^{\mu \nu }(p,p)=y^\mu y^\nu {\displaystyle \frac{pn}{py}}_{\widehat{𝐲}}.`$ (A.71)
The angular brackets indicate averaging over the spatial directions of the light-like vector $`y=(1,\widehat{𝐲})`$. The tensor $`𝒯^{\mu \nu }`$ is symmetric in $`\mu `$ and $`\nu `$ and satisfies the “Ward identity”

$`p_\mu 𝒯^{\mu \nu }(p,p)=pnn^\nu .`$ (A.72)
The self-energy tensor $`\mathrm{\Pi }^{\mu \nu }`$ is therefore also symmetric in $`\mu `$ and $`\nu `$ and satisfies

$`p_\mu \mathrm{\Pi }^{\mu \nu }\left(p\right)`$ $`=`$ $`0,`$ (A.73)
$`g_{\mu \nu }\mathrm{\Pi }^{\mu \nu }\left(p\right)`$ $`=`$ $`m_D^2.`$ (A.74)
The gluon self-energy tensor can be expressed in terms of two scalar functions, the transverse and longitudinal self-energies $`\mathrm{\Pi }_T`$ and $`\mathrm{\Pi }_L`$, defined by

$`\mathrm{\Pi }_T\left(p\right)`$ $`=`$ $`{\displaystyle \frac{1}{d1}}\left(\delta ^{ij}\widehat{p}^i\widehat{p}^j\right)\mathrm{\Pi }^{ij}\left(p\right),`$ (A.75)
$`\mathrm{\Pi }_L\left(p\right)`$ $`=`$ $`\mathrm{\Pi }^{00}\left(p\right),`$ (A.76)
where $`\widehat{𝐩}`$ is the unit vector in the direction of $`𝐩`$. In terms of these functions, the self-energy tensor is

$`\mathrm{\Pi }^{\mu \nu }\left(p\right)=\mathrm{\Pi }_T\left(p\right)T_p^{\mu \nu }{\displaystyle \frac{1}{n_p^2}}\mathrm{\Pi }_L\left(p\right)L_p^{\mu \nu },`$ (A.77)
where the tensors $`T_p`$ and $`L_p`$ are

$`T_p^{\mu \nu }`$ $`=`$ $`g^{\mu \nu }{\displaystyle \frac{p^\mu p^\nu }{p^2}}{\displaystyle \frac{n_p^\mu n_p^\nu }{n_p^2}},`$ (A.78)
$`L_p^{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{n_p^\mu n_p^\nu }{n_p^2}}.`$ (A.79)
The four-vector $`n_p^\mu `$ is

$`n_p^\mu =n^\mu {\displaystyle \frac{np}{p^2}}p^\mu `$ (A.80)
and satisfies $`pn_p=0`$ and $`n_p^2=1\left(np\right)^2/p^2`$. The equation (A.74) reduces to the identity

$`\left(d1\right)\mathrm{\Pi }_T\left(p\right)+{\displaystyle \frac{1}{n_p^2}}\mathrm{\Pi }_L\left(p\right)=m_D^2.`$ (A.81)
We can express both self-energy functions in terms of the function $`𝒯^{00}`$ defined by (A.71):

$`\mathrm{\Pi }_T\left(p\right)`$ $`=`$ $`{\displaystyle \frac{m_D^2}{\left(d1\right)n_p^2}}\left[𝒯^{00}(p,p)1+n_p^2\right],`$ (A.82)
$`\mathrm{\Pi }_L\left(p\right)`$ $`=`$ $`m_D^2\left[1𝒯^{00}(p,p)\right],`$ (A.83)
In the tensor $`𝒯^{\mu \nu }(p,p)`$ defined in (A.71), the angular brackets indicate the angular average over the unit vector $`\widehat{𝐲}`$. In almost all previous work, the angular average in (A.71) has been taken in $`d=3`$ dimensions. For consistency of higher order radiative corrections, it is essential to take the angular average in $`d=32ϵ`$ dimensions and analytically continue to $`d=3`$ only after all poles in $`ϵ`$ have been cancelled. Expressing the angular average as an integral over the cosine of an angle, the expression for the $`00`$ component of the tensor is

$`𝒯^{00}(p,p)`$ $`=`$ $`{\displaystyle \frac{w\left(ϵ\right)}{2}}{\displaystyle _1^1}𝑑c\left(1c^2\right)^ϵ{\displaystyle \frac{p_0}{p_0\left|𝐩\right|c}},`$ (A.84)
where the weight function $`w\left(ϵ\right)`$ is

$`w\left(ϵ\right)={\displaystyle \frac{\mathrm{\Gamma }\left(22ϵ\right)}{\mathrm{\Gamma }^2\left(1ϵ\right)}}2^{2ϵ}={\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3}{2}ϵ\right)}{\mathrm{\Gamma }\left(\frac{3}{2}\right)\mathrm{\Gamma }\left(1ϵ\right)}}.`$ (A.85)
The integral in (A.84) must be defined so that it is analytic at $`p_0=\mathrm{}`$. It then has a branch cut running from $`p_0=\left|𝐩\right|`$ to $`p_0=+\left|𝐩\right|`$. If we take the limit $`ϵ0`$, it reduces to

$`𝒯^{00}(p,p)`$ $`=`$ $`{\displaystyle \frac{p_0}{2\left|𝐩\right|}}\mathrm{log}{\displaystyle \frac{p_0+\left|𝐩\right|}{p_0\left|𝐩\right|}},`$ (A.86)
which is the expression that appears in the usual HTL self-energy functions.

The Feynman rule for the gluon propagator is

$`i\delta ^{ab}\mathrm{\Delta }_{\mu \nu }\left(p\right),`$ (A.87)
where the gluon propagator tensor $`\mathrm{\Delta }_{\mu \nu }`$ depends on the choice of gauge fixing. We consider two possibilities that introduce an arbitrary gauge parameter $`\xi `$: general covariant gauge and general Coulomb gauge. In both cases, the inverse propagator reduces in the limit $`\xi \mathrm{}`$ to

$`\mathrm{\Delta }_{\mathrm{}}^1\left(p\right)^{\mu \nu }`$ $`=`$ $`p^2g^{\mu \nu }+p^\mu p^\nu \mathrm{\Pi }^{\mu \nu }\left(p\right).`$ (A.88)
This can also be written

$`\mathrm{\Delta }_{\mathrm{}}^1\left(p\right)^{\mu \nu }`$ $`=`$ $`{\displaystyle \frac{1}{\mathrm{\Delta }_T\left(p\right)}}T_p^{\mu \nu }+{\displaystyle \frac{1}{n_p^2\mathrm{\Delta }_L\left(p\right)}}L_p^{\mu \nu },`$ (A.89)
where $`\mathrm{\Delta }_T`$ and $`\mathrm{\Delta }_L`$ are the transverse and longitudinal propagators:

$`\mathrm{\Delta }_T\left(p\right)`$ $`=`$ $`{\displaystyle \frac{1}{p^2\mathrm{\Pi }_T\left(p\right)}},`$ (A.90)
$`\mathrm{\Delta }_L\left(p\right)`$ $`=`$ $`{\displaystyle \frac{1}{n_p^2p^2+\mathrm{\Pi }_L\left(p\right)}}.`$ (A.91)
The inverse propagator for general $`\xi `$ is

$`\mathrm{\Delta }^1\left(p\right)^{\mu \nu }`$ $`=`$ $`\mathrm{\Delta }_{\mathrm{}}^1\left(p\right)^{\mu \nu }{\displaystyle \frac{1}{\xi }}p^\mu p^\nu \text{covariant},`$ (A.93)
$`=`$ $`\mathrm{\Delta }_{\mathrm{}}^1\left(p\right)^{\mu \nu }{\displaystyle \frac{1}{\xi }}\left(p^\mu pnn^\mu \right)\left(p^\nu pnn^\nu \right)`$
$`\text{Coulomb}.`$
The propagators obtained by inverting the tensors in (A.93) and (A.93) are

$`\mathrm{\Delta }^{\mu \nu }\left(p\right)`$ $`=`$ $`\mathrm{\Delta }_T\left(p\right)T_p^{\mu \nu }+\mathrm{\Delta }_L\left(p\right)n_p^\mu n_p^\nu \xi {\displaystyle \frac{p^\mu p^\nu }{\left(p^2\right)^2}}`$ (A.95)
$`\text{covariant},`$
$`=`$ $`\mathrm{\Delta }_T\left(p\right)T_p^{\mu \nu }+\mathrm{\Delta }_L\left(p\right)n^\mu n^\nu \xi {\displaystyle \frac{p^\mu p^\nu }{\left(n_p^2p^2\right)^2}}`$
$`\text{Coulomb}.`$
It is convenient to define the following combination of propagators:

$`\mathrm{\Delta }_X\left(p\right)`$ $`=`$ $`\mathrm{\Delta }_L\left(p\right)+{\displaystyle \frac{1}{n_p^2}}\mathrm{\Delta }_T\left(p\right).`$ (A.96)
Using (A.81), (A.90), and (A.91), it can be expressed in the alternative form

$`\mathrm{\Delta }_X\left(p\right)`$ $`=`$ $`\left[m_D^2d\mathrm{\Pi }_T\left(p\right)\right]\mathrm{\Delta }_L\left(p\right)\mathrm{\Delta }_T\left(p\right),`$ (A.97)
which shows that it vanishes in the limit $`m_D0`$. In the covariant gauge, the propagator tensor can be written

$`\mathrm{\Delta }^{\mu \nu }\left(p\right)=\left[\mathrm{\Delta }_T\left(p\right)g^{\mu \nu }+\mathrm{\Delta }_X\left(p\right)n^\mu n^\nu \right]`$
$`{\displaystyle \frac{np}{p^2}}\mathrm{\Delta }_X\left(p\right)\left(p^\mu n^\nu +n^\mu p^\nu \right)`$
$`+\left[\mathrm{\Delta }_T\left(p\right)+{\displaystyle \frac{\left(np\right)^2}{p^2}}\mathrm{\Delta }_X\left(p\right){\displaystyle \frac{\xi }{p^2}}\right]{\displaystyle \frac{p^\mu p^\nu }{p^2}}.`$ (A.98)
This decomposition of the propagator into three terms has proved to be particularly convenient for explicit calculations. For example, the first term satisfies the identity

$`\left[\mathrm{\Delta }_T\left(p\right)g_{\mu \nu }+\mathrm{\Delta }_X\left(p\right)n_\mu n_\nu \right]\mathrm{\Delta }_{\mathrm{}}^1\left(p\right)^{\nu \lambda }=`$
$`g_{\mu }^{}{}_{}{}^{\lambda }{\displaystyle \frac{p_\mu p^\lambda }{p^2}}+{\displaystyle \frac{np}{n_p^2p^2}}{\displaystyle \frac{\mathrm{\Delta }_X\left(p\right)}{\mathrm{\Delta }_L\left(p\right)}}p_\mu n_p^\lambda .`$ (A.99)
### A.2 Quark Self-energy

The HTL self-energy of a quark with momentum $`p`$ is given by

$`\mathrm{\Sigma }\left(P\right)=m_q^2/𝒯\left(p\right),`$ (A.100)
where

$`𝒯^\mu \left(p\right)={\displaystyle \frac{y^\mu }{py}}_{\widehat{𝐲}}.`$ (A.101)
Expressing the angular average as an integral over the cosine of an angle, the expression is

$`𝒯^\mu \left(p\right)={\displaystyle \frac{w\left(ϵ\right)}{2}}{\displaystyle _1^1}𝑑c\left(1c^2\right)^ϵ{\displaystyle \frac{y^\mu }{p_0\left|𝐩\right|c}},`$ (A.102)
The integral in (A.102) must be defined so that it is analytic at $`p_0=\mathrm{}`$. It then has a branch cut running from $`p_0=\left|𝐩\right|`$ to $`p_0=+\left|𝐩\right|`$. In three dimensions, this reduces to

$`\mathrm{\Sigma }\left(P\right)`$ $`=`$ $`{\displaystyle \frac{m_q^2}{2\left|𝐩\right|}}\gamma _0\mathrm{log}{\displaystyle \frac{p_0+\left|𝐩\right|}{p_0\left|𝐩\right|}}`$ (A.103)
$`+{\displaystyle \frac{m_q^2}{\left|𝐩\right|}}\gamma \widehat{𝐩}\left(1{\displaystyle \frac{p_0}{2\left|𝐩\right|}}\mathrm{log}{\displaystyle \frac{p_0+\left|𝐩\right|}{p_0\left|𝐩\right|}}\right).`$
### A.3 Quark Propagator

The Feynman rule for the quark propagator is

$`i\delta ^{ab}S\left(p\right).`$ (A.104)
The quark propagator can be written as

$`S\left(p\right)={\displaystyle \frac{1}{/p\mathrm{\Sigma }\left(p\right)}},`$ (A.105)
where the quark self-energy is given by (A.100). The inverse quark propagator can be written as

$`S^1\left(p\right)=/p\mathrm{\Sigma }\left(p\right).`$ (A.106)
This can be written as

$`S^1\left(p\right)=/𝒜\left(p\right),`$ (A.107)
where we have organized $`A_0\left(p\right)`$ and $`A_S\left(p\right)`$ into:

$`A_\mu \left(p\right)=(A_0\left(p\right),A_S\left(p\right)\widehat{𝐩}).`$ (A.108)
The functions $`A_0\left(p\right)`$ and $`A_S\left(p\right)`$ are defined

$`A_0\left(p\right)`$ $`=`$ $`p_0{\displaystyle \frac{m_q^2}{p_0}}𝒯_p,`$ (A.109)
$`A_S\left(p\right)`$ $`=`$ $`\left|𝐩\right|+{\displaystyle \frac{m_q^2}{\left|𝐩\right|}}\left[1𝒯_p\right].`$ (A.110)
### A.4 Quark-gluon vertex

The quark-gluon vertex with outgoing gluon momentum $`p`$, incoming fermion momentum $`q`$, and outgoing quark momentum $`r`$, Lorentz index $`\mu `$ and color index $`a`$ is

$`\mathrm{\Gamma }_a^\mu (p,q,r)=gt_a\left(\gamma ^\mu m_q^2\stackrel{~}{𝒯}^\mu (p,q,r)\right).`$ (A.111)
The tensor in the HTL correction term is only defined for $`pq+r=0`$:

$`\stackrel{~}{𝒯}^\mu (p,q,r)=y^\mu \left({\displaystyle \frac{y/}{qyry}}\right)_{\widehat{𝐲}}.`$ (A.112)
This tensor is even under the permutation of $`q`$ and $`r`$. It satisfies the “Ward identity”

$`p_\mu \stackrel{~}{𝒯}^\mu (p,q,r)=\stackrel{~}{𝒯}^\mu \left(q\right)\stackrel{~}{𝒯}^\mu \left(r\right).`$ (A.113)
The quark-gluon vertex therefore satisfies the Ward identity

$`p_\mu \mathrm{\Gamma }^\mu (p,q,r)=S^1\left(q\right)S^1\left(r\right).`$ (A.114)
### A.5 Quark-gluon four-vertex

We define the quark-gluon four-point vertex with outgoing gluon momenta $`p`$ and $`q`$, incoming fermion momentum $`r`$, and outgoing fermion momentum $`s`$. Generally this vertex has both adjoint and fundamental indices, however, for this calculation we will only need the quark-gluon four-point vertex traced over the adjoint color indices. In this case

$`\delta ^{ab}\mathrm{\Gamma }_{abij}^{\mu \nu }(p,q,r,s)`$ $`=`$ $`g^2m_q^2c_F\delta _{ij}\stackrel{~}{𝒯}^{\mu \nu }(p,q,r,s)`$ (A.115)
$``$ $`g^2c_F\delta _{ij}\mathrm{\Gamma }^{\mu \nu },`$
where $`c_F=\left(N_c^21\right)/\left(2N_c\right)`$. There is no tree-level term. The tensor in the HTL correction term is only defined for $`p+qr+s=0`$

$`\stackrel{~}{𝒯}^{\mu \nu }(p,q,r,s)`$ $`=`$ $`y^\mu y^\nu ({\displaystyle \frac{1}{ry}}+{\displaystyle \frac{1}{sy}})`$ (A.116)
$`\times {\displaystyle \frac{y/}{\left[\left(rp\right)y\right]\left[\left(s+p\right)y\right]}}.`$
This tensor is symmetric in $`\mu `$ and $`\nu `$ and is traceless. It satisfies the Ward identity:

$`p_\mu \mathrm{\Gamma }^{\mu \nu }(p,q,r,s)=\mathrm{\Gamma }^\nu (q,rp,s)\mathrm{\Gamma }^\nu (q,r,s+p).`$ (A.117)
### A.6 HTL Quark Counterterm

The Feynman rule for the insertion of an HTL quark counterterm into a quark propagator is

$`i\delta ^{ab}\mathrm{\Sigma }\left(p\right),`$ (A.118)
where $`\mathrm{\Sigma }\left(p\right)`$ is the HTL quark self-energy given in (A.108).

### A.7 Imaginary-time formalism

In the imaginary-time formalism, Minkoswski energies have discrete imaginary values $`p_0=i\left(2\pi nT\right)`$ and integrals over Minkowski space are replaced by sum-integrals over Euclidean vectors $`(2\pi nT,𝐩)`$. We will use the notation $`P=(P_0,𝐩)`$ for Euclidean momenta. The magnitude of the spatial momentum will be denoted $`p=\left|𝐩\right|`$, and should not be confused with a Minkowski vector. The inner product of two Euclidean vectors is $`PQ=P_0Q_0+𝐩𝐪`$. The vector that specifies the thermal rest frame remains $`n=(1,\mathrm{𝟎})`$.

The Feynman rules for Minkowski space given above can be easily adapted to Euclidean space. The Euclidean tensor in a given Feynman rule is obtained from the corresponding Minkowski tensor with raised indices by replacing each Minkowski energy $`p_0`$ by $`iP_0`$, where $`P_0`$ is the corresponding Euclidean energy, and multipying by $`i`$ for every $`0`$ index. This prescription transforms $`p=(p_0,𝐩)`$ into $`P=(P_0,𝐩)`$, $`g^{\mu \nu }`$ into $`\delta ^{\mu \nu }`$, and $`pq`$ into $`PQ`$. The effect on the HTL tensors defined in (A.71), (A.112), and (A.116) is equivalent to substituting $`pnPN`$ where $`N=(i,\mathrm{𝟎})`$, $`pyPY`$ where $`Y=(i,\widehat{𝐲})`$, and $`y^\mu Y^\mu `$. For example, the Euclidean tensor corresponding to (A.71) is

$`𝒯^{\mu \nu }(P,P)=Y^\mu Y^\nu {\displaystyle \frac{PN}{PY}}.`$ (A.119)
The average is taken over the directions of the unit vector $`\widehat{𝐲}`$.

Alternatively, one can calculate a diagram by using the Feynman rules for Minkowski momenta, reducing the expressions for diagrams to scalars, and then make the appropriate substitutions, such as $`p^2P^2`$, $`pqPQ`$, and $`npinP`$. For example, the propagator functions (A.90) and (A.91) become

$`\mathrm{\Delta }_T\left(P\right)`$ $`=`$ $`{\displaystyle \frac{1}{P^2+\mathrm{\Pi }_T\left(P\right)}},`$ (A.120)
$`\mathrm{\Delta }_L\left(P\right)`$ $`=`$ $`{\displaystyle \frac{1}{p^2+\mathrm{\Pi }_L\left(P\right)}}.`$ (A.121)
The expressions for the HTL self-energy functions $`\mathrm{\Pi }_T\left(P\right)`$ and $`\mathrm{\Pi }_L\left(P\right)`$ are given by (A.82) and (A.83) with $`n_p^2`$ replaced by $`n_P^2=p^2/P^2`$ and $`𝒯^{00}(p,p)`$ replaced by

$`𝒯_P`$ $`=`$ $`{\displaystyle \frac{w\left(ϵ\right)}{2}}{\displaystyle _1^1}𝑑c\left(1c^2\right)^ϵ{\displaystyle \frac{iP_0}{iP_0pc}}.`$ (A.122)
Note that this function differs by a sign from the 00 component of the Euclidean tensor corresponding to (A.71):

$`𝒯^{00}(P,P)=𝒯^{00}(p,p)|_{p_0iP_0}=𝒯_P.`$ (A.123)
A more convenient form for calculating sum-integrals that involve the function $`𝒯_P`$ is

$`𝒯_P`$ $`=`$ $`{\displaystyle \frac{P_0^2}{P_0^2+p^2c^2}}_c,`$ (A.124)
where the angular brackets represent an average over $`c`$ defined by

$$f\left(c\right)_cw\left(ϵ\right)_0^1𝑑c\left(1c^2\right)^ϵf\left(c\right)$$
(A.125)
and $`w\left(ϵ\right)`$ is given in (A.85).

## Appendix B Sum-integrals

In the imaginary-time formalism for thermal field theory, the 4-momentum $`P=(P_0,𝐩)`$is Euclidean with $`P^2=P_0^2+𝐩^2`$. The Euclidean energy $`p_0`$ has discrete values: $`P_0=2n\pi T`$ for bosons and $`P_0=\left(2n+1\right)\pi T`$ for fermions, where $`n`$ is an integer. Loop diagrams involve sums over $`P_0`$ and integrals over $`𝐩`$. With dimensional regularization, the integral is generalized to $`d=32ϵ`$ spatial dimensions. We define the dimensionally regularized sum-integral by

$`{\displaystyle _P}`$ $``$ $`\left({\displaystyle \frac{e^\gamma \mu ^2}{4\pi }}\right)^ϵT{\displaystyle \underset{P_0=2n\pi T}{}}{\displaystyle \frac{d^{32ϵ}p}{\left(2\pi \right)^{32ϵ}}},`$ (B.1)
$`{\displaystyle _{\left\{P\right\}}}`$ $``$ $`\left({\displaystyle \frac{e^\gamma \mu ^2}{4\pi }}\right)^ϵT{\displaystyle \underset{P_0=\left(2n+1\right)\pi T}{}}{\displaystyle \frac{d^{32ϵ}p}{\left(2\pi \right)^{32ϵ}}},`$ (B.2)
where $`32ϵ`$ is the dimension of space and $`\mu `$ is an arbitrary momentum scale. The factor $`\left(e^\gamma /4\pi \right)^ϵ`$ is introduced so that, after minimal subtraction of the poles in $`ϵ`$ due to ultraviolet divergences, $`\mu `$ coincides with the renormalization scale of the $`\overline{\mathrm{MS}}`$ renormalization scheme.

### B.1 One-loop sum-integrals

The simple one-loop sum-integrals required in our calculations can be derived from the formulas

$`{\displaystyle _P}{\displaystyle \frac{p^{2m}}{\left(P^2\right)^n}}`$ $`=`$ $`\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}{\displaystyle \frac{2\mathrm{\Gamma }\left(\frac{3}{2}+mϵ\right)\mathrm{\Gamma }\left(n\frac{3}{2}m+ϵ\right)}{\mathrm{\Gamma }\left(n\right)\mathrm{\Gamma }\left(22ϵ\right)}}`$ (B.3)
$`\times \mathrm{\Gamma }\left(1ϵ\right)\zeta \left(2n2m3+2ϵ\right)e^{ϵ\gamma }`$
$`\times T^{4+2m2n}\left(2\pi \right)^{1+2m2n},`$
$`{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{p^{2m}}{\left(P^2\right)^n}}`$ $`=`$ $`\left(2^{2n2md}1\right){\displaystyle _P}{\displaystyle \frac{p^{2m}}{\left(P^2\right)^n}}.`$ (B.4)
The specific bosonic one-loop sum-integrals needed are

$`{\displaystyle _P}{\displaystyle \frac{1}{P^2}}`$ $`=`$ $`{\displaystyle \frac{T^2}{12}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}[1+(2+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}})ϵ`$ (B.5)
$`+(4+{\displaystyle \frac{\pi ^2}{4}}+4{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}+2{\displaystyle \frac{\zeta ^{\prime \prime }\left(1\right)}{\zeta \left(1\right)}})ϵ^2],`$
$`{\displaystyle _P}{\displaystyle \frac{p^2}{\left(P^2\right)^2}}`$ $`=`$ $`{\displaystyle \frac{1}{8}}T^2,`$ (B.6)
$`{\displaystyle _P}{\displaystyle \frac{1}{\left(P^2\right)^2}}`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}`$ (B.7)
$`\times \left[{\displaystyle \frac{1}{ϵ}}+2\gamma +\left({\displaystyle \frac{\pi ^2}{4}}4\gamma _1\right)ϵ\right],`$
$`{\displaystyle _P}{\displaystyle \frac{1}{p^2P^2}}`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}2[{\displaystyle \frac{1}{ϵ}}+2\gamma +2`$ (B.8)
$`+(4+4\gamma +{\displaystyle \frac{\pi ^2}{4}}4\gamma _1)ϵ].`$
The specific fermionic one-loop sum-integrals needed are

$`{\displaystyle _{\left\{P\right\}}}\mathrm{log}P^2`$ $`=`$ $`{\displaystyle \frac{7\pi ^2}{360}}T^4,`$ (B.9)
$`{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{1}{P^2}}`$ $`=`$ $`{\displaystyle \frac{T^2}{24}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}`$ (B.10)
$`\times [1+(22\mathrm{log}2+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}})ϵ`$
$`+(4+{\displaystyle \frac{\pi ^2}{4}}4\mathrm{log}22\mathrm{log}^22`$
$`+4(1\mathrm{log}2){\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}+2{\displaystyle \frac{\zeta ^{\prime \prime }\left(1\right)}{\zeta \left(1\right)}})ϵ^2],`$
$`{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{1}{\left(P^2\right)^2}}`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}\left[{\displaystyle \frac{1}{ϵ}}+2\gamma +4\mathrm{log}2\right],`$ (B.11)
$`{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{p^2}{\left(P^2\right)^2}}`$ $`=`$ $`{\displaystyle \frac{T^2}{16}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}`$ (B.12)
$`\times \left[1+\left({\displaystyle \frac{4}{3}}2\mathrm{log}2+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right)ϵ\right],`$
$`{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{p^2}{\left(P^2\right)^3}}`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}`$ (B.13)
$`\times {\displaystyle \frac{3}{4}}\left[{\displaystyle \frac{1}{ϵ}}+2\gamma {\displaystyle \frac{2}{3}}+4\mathrm{log}2\right],`$
$`{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{p^4}{\left(P^2\right)^3}}`$ $`=`$ $`{\displaystyle \frac{5T^2}{64}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}`$ (B.14)
$`\times \left[1+\left({\displaystyle \frac{14}{15}}2\mathrm{log}2+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right)ϵ\right],`$
$`{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{p^4}{\left(P^2\right)^4}}`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}`$ (B.15)
$`\times {\displaystyle \frac{5}{8}}\left[{\displaystyle \frac{1}{ϵ}}+2\gamma {\displaystyle \frac{16}{15}}+4\mathrm{log}2\right],`$
$`{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{1}{p^2P^2}}`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}2[{\displaystyle \frac{1}{ϵ}}+2+2\gamma +4\mathrm{log}2`$ (B.16)
$`+(4+8\mathrm{log}2+4\mathrm{log}^22+4\gamma (1+2\mathrm{log}2)`$
$`+{\displaystyle \frac{\pi ^2}{4}}4\gamma _1)ϵ].`$
The errors are all of one order higher in $`ϵ`$ than the smallest term shown. The number $`\gamma _1`$ is the first Stieltjes gamma constant defined by the equation

$$\zeta \left(1+z\right)=\frac{1}{z}+\gamma \gamma _1z+O\left(z^2\right).$$
(B.17)
### B.2 One-loop HTL sum-integrals

We also need some more difficult one-loop sum-integrals that involve the HTL function defined in (A.102).

The specific bosonic sum-integrals needed are

$`{\displaystyle _P}{\displaystyle \frac{1}{p^4}}𝒯_P`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}`$ (B.18)
$`\times \left(1\right)\left[{\displaystyle \frac{1}{ϵ}}+2\gamma +2\mathrm{log}2\right],`$
$`{\displaystyle _P}{\displaystyle \frac{1}{p^2P^2}}𝒯_P`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}`$ (B.19)
$`\times \left[2\mathrm{log}2\left({\displaystyle \frac{1}{ϵ}}+2\gamma \right)+2\mathrm{log}^22+{\displaystyle \frac{\pi ^2}{3}}\right],`$
$`{\displaystyle _P}{\displaystyle \frac{1}{\left(P^2\right)^2}}𝒯_P`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}{\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{1}{ϵ}}+2\gamma +1\right].`$ (B.20)
The specific fermionic sum-integrals needed are

$`{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{1}{\left(P^2\right)^2}}𝒯_P`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}`$ (B.21)
$`\times {\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{1}{ϵ}}+2\gamma +1+4\mathrm{log}2\right],`$
$`{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{1}{p^2P^2}}𝒯_P`$ $`=`$ $`{\displaystyle \frac{2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}`$ (B.22)
$`\times \left[\mathrm{log}2\left({\displaystyle \frac{1}{ϵ}}+2\gamma \right)+5\mathrm{log}^22+{\displaystyle \frac{\pi ^2}{6}}\right],`$
$`{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{1}{P^2P_0^2}}𝒯_P`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}`$ (B.23)
$`\times [{\displaystyle \frac{1}{ϵ^2}}+2(\gamma +2\mathrm{log}2){\displaystyle \frac{1}{ϵ}}+{\displaystyle \frac{\pi ^2}{4}}`$
$`+4\mathrm{log}^22+8\gamma \mathrm{log}24\gamma _1],`$
$`{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{1}{p^2P_0^2}}\left(𝒯_P\right)^2`$ $`=`$ $`{\displaystyle \frac{4}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}`$ (B.24)
$`\times \left[\mathrm{log}2\left({\displaystyle \frac{1}{ϵ}}+2\gamma \right)+5\mathrm{log}^22\right],`$
$`{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{1}{P^2}}{\displaystyle \frac{1}{\left(QY\right)^2}}_c`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{2ϵ}`$ (B.25)
$`\times \left(1\right)\left[{\displaystyle \frac{1}{ϵ}}1+2\gamma +4\mathrm{log}2\right].`$
The errors are all of order $`ϵ`$

It is straightforward to calculate the sum-integrals (B.21)–(B.24) using the representation (A.124) of the function $`𝒯_P`$. For example, the sum-integral (B.18) can be written

$`{\displaystyle _P}{\displaystyle \frac{1}{p^4}}𝒯_P={\displaystyle _P}{\displaystyle \frac{1}{p^4}}{\displaystyle \frac{P_0^2}{P_0^2+p^2c^2}}_c,`$ (B.26)
where the angular brackets denote an average over $`c`$ as defined in (A.125).

$`{\displaystyle _P}{\displaystyle \frac{1}{p^4}}𝒯_P={\displaystyle _P}{\displaystyle \frac{1}{p^4}}\left[1{\displaystyle \frac{p^2c^2}{P_0^2+p^2c^2}}_c\right].`$ (B.27)
The first term in the square brackets vanishes with dimensional regularization, while after rescaling the momentum by $`𝐩𝐩/c`$, the second term reads

$`{\displaystyle _P}{\displaystyle \frac{1}{p^4}}𝒯_P=c^{1+2ϵ}_c{\displaystyle _P}{\displaystyle \frac{1}{p^2P^2}}.`$ (B.28)
Evaluating the average over $`c`$, using the expression (B.8) for the sum-integral, and expanding in powers of $`ϵ`$, we obtain the result (B.18). Following the same strategy, all the sum-integrals (B.21)–(B.24) can be reduced to linear combinations of simple sum-integrals with coefficients that are averages over $`c`$. The only difficult integral is the double average over $`c`$ that arises from (B.24):

$`{\displaystyle \frac{c_1^{1+2ϵ}c_2^{1+2ϵ}}{c_1^2c_2^2}}_{c_1,c_2}=2\mathrm{log}2+2\left(\mathrm{log}^222\mathrm{log}2\right)ϵ.`$
### B.3 Simple two-loop sum-integrals

The simple two-loop sum-integrals that are needed are

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{1}{P^2Q^2R^2}}`$ $`=`$ $`0,`$ (B.30)
$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{1}{P^2Q^2r^2}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$ (B.31)
$`\times \left({\displaystyle \frac{1}{6}}\right)\left[{\displaystyle \frac{1}{ϵ}}+42\mathrm{log}2+4{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right],`$
$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{q^2}{P^2Q^2r^4}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$ (B.32)
$`\times \left({\displaystyle \frac{1}{12}}\right)\left[{\displaystyle \frac{1}{ϵ}}+{\displaystyle \frac{11}{3}}+2\gamma 2\mathrm{log}2+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right],`$
$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{q^2}{P^2Q^2r^2R^2}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$ (B.33)
$`\times \left({\displaystyle \frac{1}{72}}\right)\left[{\displaystyle \frac{1}{ϵ}}7.002\right],`$
$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{PQ}{P^2Q^2r^4}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$ (B.34)
$`\times \left({\displaystyle \frac{1}{36}}\right)\left[16\gamma +6{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right],`$
$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{p^2}{q^2P^2Q^2R^2}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$ (B.35)
$`\times \left({\displaystyle \frac{5}{72}}\right)\left[{\displaystyle \frac{1}{ϵ}}+9.5667\right],`$
$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{r^2}{q^2P^2Q^2R^2}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$ (B.36)
$`\times \left({\displaystyle \frac{1}{18}}\right)\left[{\displaystyle \frac{1}{ϵ}}+8.1420\right],`$
where $`R=\left(P+Q\right)`$ and $`r=\left|𝐩+𝐪\right|`$. The corrections are all of order $`ϵ`$. To motivate the integration formula we will use to evaluate the two-loop sum-integrals, we first present the analogous integration formula for one-loop sum-integrals. In a one-loop sum-integral, the sum over $`P_0`$ can be replaced by a contour integral in $`p_0=iP_0`$:

$`{\displaystyle _P}F\left(P\right)`$ $`=`$ $`\underset{\eta 0^+}{lim}{\displaystyle \frac{dp_0}{2\pi i}_𝐩\left[F(ip_0,𝐩)F(0,𝐩)\right]}`$ (B.37)
$`\times e^{\eta p_0}n\left(p_0\right),`$
where $`n\left(p_0\right)=1/\left(e^{\beta p_0}1\right)`$ is the Bose-Einstein thermal distribution and the contour runs from $`\mathrm{}`$ to $`+\mathrm{}`$ above the real axis and from $`+\mathrm{}`$ to $`\mathrm{}`$ below the real axis. This formula can be expressed in a more convenient form by collapsing the contour onto the real axis and separating out those terms with the exponential convergence factor $`n\left(\left|p_0\right|\right)`$. The remaining terms run along contours from $`\mathrm{}\pm i\epsilon `$ to 0 and have the convergence factor $`e^{\eta p_0}`$. This allows the contours to be deformed so that they run from 0 to $`\pm i\mathrm{}`$ along the imaginary $`p_0`$ axis, which corresponds to real values of $`P_0=ip_0`$. Assuming that $`F(ip_0,𝐩)`$ is a real function of $`p_0`$, i.e. that it satisfies $`F(ip_0^{},𝐩)=F(ip_0,𝐩)^{}`$, the resulting formula for the sum-integral is

$`{\displaystyle _P}F\left(P\right)`$ $`=`$ $`{\displaystyle _P}F\left(P\right)`$ (B.38)
$`+{\displaystyle _p}ϵ\left(p_0\right)n\left(\left|p_0\right|\right)\mathrm{\hspace{0.17em}2}\mathrm{Im}F(ip_0+\epsilon ,𝐩),`$
where $`ϵ\left(p_0\right)`$ is the sign of $`p_0`$. The first integral on the right side is over the $`\left(d+1\right)`$-dimensional Euclidean vector $`P=(P_0,𝐩)`$ and the second is over the $`\left(d+1\right)`$-dimensional Minkowskian vector $`p=(p_0,𝐩)`$.

The two-loop sum-integrals can be evaluated by using a generalization of the one-loop formula (B.38):

$`{\displaystyle _{\left\{PQ\right\}}}F\left(P\right)G\left(Q\right)H\left(R\right)`$ $`=`$ $`{\displaystyle _{PQ}}F\left(P\right)G\left(Q\right)H\left(R\right){\displaystyle _p}ϵ\left(p_0\right)n_F\left(\left|p_0\right|\right)\mathrm{\hspace{0.17em}2}\mathrm{Im}F(ip_0+\epsilon ,𝐩)\mathrm{Re}{\displaystyle _Q}G\left(Q\right)H\left(R\right)|_{P_0=ip_0+\epsilon }`$ (B.39)
$`{\displaystyle _p}ϵ\left(p_0\right)n_F\left(\left|p_0\right|\right)\mathrm{\hspace{0.17em}2}\mathrm{Im}G(ip_0+\epsilon ,𝐩)\mathrm{Re}{\displaystyle _Q}H\left(Q\right)F\left(R\right)|_{P_0=ip_0+\epsilon }`$
$`+{\displaystyle _p}ϵ\left(p_0\right)n_B\left(\left|p_0\right|\right)\mathrm{\hspace{0.17em}2}\mathrm{Im}H(ip_0+\epsilon ,𝐩)\mathrm{Re}{\displaystyle _Q}F\left(Q\right)G\left(R\right)|_{P_0=ip_0+\epsilon }`$
$`+{\displaystyle _p}ϵ\left(p_0\right)n_F\left(\left|p_0\right|\right)\mathrm{\hspace{0.17em}2}\mathrm{Im}F(ip_0+\epsilon ,𝐩){\displaystyle _q}ϵ\left(q_0\right)n_F\left(\left|q_0\right|\right)\mathrm{\hspace{0.17em}2}\mathrm{Im}G(iq_0+\epsilon ,𝐪)\mathrm{Re}H\left(R\right)|_{R_0=i\left(p_0+q_0\right)+\epsilon }`$
$`{\displaystyle _p}ϵ\left(p_0\right)n_F\left(\left|p_0\right|\right)\mathrm{\hspace{0.17em}2}\mathrm{Im}G(ip_0+\epsilon ,𝐩){\displaystyle _q}ϵ\left(q_0\right)n_B\left(\left|q_0\right|\right)\mathrm{\hspace{0.17em}2}\mathrm{Im}H(iq_0+\epsilon ,𝐪)\mathrm{Re}F\left(R\right)|_{R_0=i\left(p_0+q_0\right)+\epsilon }`$
$`{\displaystyle _p}ϵ\left(p_0\right)n_B\left(\left|p_0\right|\right)\mathrm{\hspace{0.17em}2}\mathrm{Im}H(ip_0+\epsilon ,𝐩){\displaystyle _q}ϵ\left(q_0\right)n_F\left(\left|q_0\right|\right)\mathrm{\hspace{0.17em}2}\mathrm{Im}F(iq_0+\epsilon ,𝐪)\mathrm{Re}G\left(R\right)|_{R_0=i\left(p_0+q_0\right)+\epsilon }.`$
This formula can be derived in 3 steps. First, express the sum over $`P_0`$ as the sum of two contour integrals over $`p_0`$, one that encloses the real axis $`\mathrm{Im}p_0=0`$ and another that encloses the line $`\mathrm{Im}p_0=\mathrm{Im}q_0`$. Second, express the the sum over $`q_0`$ as a contour integral that encloses the real-$`q_0`$ axis. The resulting terms can be combined into the expression (B.39). The integrals of the imaginary parts that enter into our calculation can be reduced to

$`{\displaystyle _p}ϵ\left(p_0\right)n\left(\left|p_0\right|\right)\mathrm{\hspace{0.17em}2}\mathrm{Im}{\displaystyle \frac{1}{P^2}}|_{P_0=ip_0+\epsilon }f(ip_0+\epsilon ,𝐩)`$
$`={\displaystyle _𝐩}{\displaystyle \frac{n\left(p\right)}{p}}{\displaystyle \frac{1}{2}}{\displaystyle \underset{\pm }{}}f(\pm ip+\epsilon ,𝐩),`$ (B.40)
$`{\displaystyle _p}ϵ\left(p_0\right)n\left(\left|p_0\right|\right)\mathrm{\hspace{0.17em}2}\mathrm{Im}𝒯_P|_{P_0=ip_0+\epsilon }f(ip_0+\epsilon ,𝐩)`$
$`={\displaystyle _𝐩}pn\left(p\right){\displaystyle \frac{1}{2}}{\displaystyle \underset{\pm }{}}c^{3+2ϵ}f(\pm ip+\epsilon ,𝐩/c)_c.`$ (B.41)
The latter equation is obtained by inserting the expression (A.124) for $`𝒯_P`$, using (B.40), and then making the change of variable $`𝐩𝐩/c`$ to put the thermal integral into a standard form.

As a simple illustration, we apply the formula (B.39) to the sum-integral (B.31). The nonvanishing terms are

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{1}{P^2Q^2r^2}}`$ $`=`$ $`2{\displaystyle _p}n_F\left(\left|p_0\right|\right)\mathrm{\hspace{0.17em}2}\pi \delta \left(p_0^2p^2\right){\displaystyle _Q}{\displaystyle \frac{1}{Q^2r^2}}`$ (B.42)
$`+{\displaystyle _p}n_F\left(\left|p_0\right|\right)\mathrm{\hspace{0.17em}2}\pi \delta \left(p_0^2p^2\right)`$
$`\times {\displaystyle _q}n_F\left(\right|q_0\left|\right)\mathrm{\hspace{0.17em}2}\pi \delta (q_0^2q^2){\displaystyle \frac{1}{r^2}}.`$
The delta functions can be used to evaluate the integrals over $`p_0`$ and $`q_0`$. The integral over $`Q`$ is given in (C.177) up to corrections of order $`ϵ`$. This reduces the sum-integral to

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{1}{P^2Q^2r^2}}`$ $`=`$ $`{\displaystyle \frac{4}{\left(4\pi \right)^2}}\left[{\displaystyle \frac{1}{ϵ}}+42\mathrm{log}2\right]\mu ^{2ϵ}`$ (B.43)
$`\times {\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}p^{2ϵ}+{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_F\left(q\right)}{pq}}{\displaystyle \frac{1}{r^2}}.`$
The momentum integrals are evaluated in (C.68) and (C.69). Keeping all terms that contribute through order $`ϵ^0`$, we get the result (B.31). The sum-integral (B.32) can be evaluated in the same way:

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{q^2}{P^2Q^2r^4}}`$ $`=`$ $`{\displaystyle \frac{2}{\left(4\pi \right)^2}}\left[{\displaystyle \frac{1}{ϵ}}2\mathrm{log}2\right]\mu ^{2ϵ}{\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}p^{2ϵ}`$ (B.44)
$`+{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_F\left(q\right)}{pq}}{\displaystyle \frac{q^2}{r^4}}.`$
The sum-integral (B.34) can be reduced to a linear combination of (B.31) and (B.32) by expressing the numerator in the form $`PQ=P_0Q_0+\left(r^2p^2q^2\right)/2`$ and noting that the $`P_0Q_0`$ term vanishes upon summing over $`P_0`$ or $`Q_0`$.

The sum-integral (B.33) is a little more difficult. After applying the formula (B.39) and using the delta functions to integrate over $`p_0`$, $`q_0`$, and $`r_0`$, it can be reduced to

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{q^2}{P^2Q^2r^2R^2}}={\displaystyle _𝐩}{\displaystyle \frac{n_B\left(p\right)}{p}}{\displaystyle _Q}{\displaystyle \frac{q^2}{p^2Q^2R^2}}|_{P=ip}`$
$`{\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle _Q}{\displaystyle \frac{1}{Q^2R^2}}\left({\displaystyle \frac{q^2}{r^2}}+{\displaystyle \frac{p^2}{q^2}}\right)|_{P=ip}`$
$`+{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_F\left(q\right)}{pq}}{\displaystyle \frac{p^2}{r^2}}{\displaystyle \frac{r^2p^2q^2}{\mathrm{\Delta }(p,q,r)}}`$
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_B\left(q\right)}{pq}}\left({\displaystyle \frac{p^2}{q^2}}+{\displaystyle \frac{r^2}{q^2}}\right){\displaystyle \frac{r^2p^2q^2}{\mathrm{\Delta }(p,q,r)}},`$ (B.45)
where $`\mathrm{\Delta }(p,q,r)`$ is the triangle function that is negative when $`p`$, $`q`$, and $`r`$ are the lengths of 3 sides of a triangle:

$$\mathrm{\Delta }(p,q,r)=p^4+q^4+r^42\left(p^2q^2+q^2r^2+r^2p^2\right).$$
(B.46)
After using (C.183)–(C.3) to integrate over $`Q`$, the first term on the right side of (B.45) is evaluated using (C.68). The 2-loop thermal integrals on the right side of (B.45) are given in (C.73)–(C.76). Adding together all the terms, we get the final result (B.33). The sum-integrals (B.35) and (B.36) are evaluated in a similar manner.

### B.4 Two-loop HTL sum-integrals

We also need some more difficult two-loop sum-integrals that involve the functions $`𝒯_P`$ defined in (A.102)

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{1}{P^2Q^2r^2}}𝒯_R`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{1}{48}}\right)\left[{\displaystyle \frac{1}{ϵ^2}}+\left(2+12\mathrm{log}2+4{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right){\displaystyle \frac{1}{ϵ}}+136.362\right],`$ (B.47)
$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{q^2}{P^2Q^2r^4}}𝒯_R`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{1}{576}}\right)\left[{\displaystyle \frac{1}{ϵ^2}}+\left({\displaystyle \frac{26}{3}}+52\mathrm{log}2+4{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right){\displaystyle \frac{1}{ϵ}}+446.438\right],`$ (B.48)
$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{PQ}{P^2Q^2r^4}}𝒯_R`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{1}{96}}\right)\left[{\displaystyle \frac{1}{ϵ^2}}+\left(4\mathrm{log}2+4{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right){\displaystyle \frac{1}{ϵ}}+69.174\right],`$ (B.49)
$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{r^2p^2}{P^2q^2Q_0^2R^2}}𝒯_Q`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}{\displaystyle \frac{1}{8}}\left[{\displaystyle \frac{1}{ϵ^2}}+\left(2+2\gamma +{\displaystyle \frac{10}{3}}\mathrm{log}2+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right){\displaystyle \frac{1}{ϵ}}+46.8757\right].`$ (B.50)
The errors are all of order $`ϵ`$. To calculate the sum-integral (B.47), we begin by using the representation (A.124) of the function $`T_R`$:

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{1}{P^2Q^2r^2}}𝒯_R`$ $`=`$ $`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{1}{P^2Q^2r^2}}`$ (B.51)
$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{1}{P^2Q^2}}{\displaystyle \frac{c^2}{R_0^2+r^2c^2}}_c.`$
The first sum-integral on the right hand side is given by (B.31). To evaluate the second sum-integral, we apply the sum-integral formula (B.39):

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{1}{P^2Q^2\left(R_0^2+r^2c^2\right)}}`$ (B.52)
$`={\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}2\mathrm{R}\mathrm{e}{\displaystyle _Q}{\displaystyle \frac{1}{Q^2\left(R_0^2+r^2c^2\right)}}|_{P_0=ip+\epsilon }`$
$`+c^{3+2ϵ}{\displaystyle _𝐩}{\displaystyle \frac{n_B\left(p\right)}{p}}{\displaystyle _Q}{\displaystyle \frac{1}{Q^2R^2}}|_{P(ip,𝐩/c)}`$
$`+{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_F\left(q\right)}{pq}}\mathrm{Re}{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}`$
$`2c^{3+2ϵ}{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_B\left(q\right)}{pq}}\mathrm{Re}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}},`$
where $`r_c=\left|𝐩+𝐪/c\right|`$. In the terms on the right side with a single thermal integral, the appropriate averages over $`c`$ of the integrals over $`Q`$ are given in (C.181) and (C.189).

The subsequent integrals over $`𝐩`$ are special cases of (C.68) and (C.69):

$`{\displaystyle _𝐩}n_B\left(p\right)p^{12ϵ}`$ $`=`$ $`2^{8ϵ}{\displaystyle \frac{\left(1\right)_{4ϵ}\left(\frac{1}{2}\right)_{2ϵ}}{\left(1\right)_{2ϵ}\left(\frac{3}{2}\right)_ϵ}}{\displaystyle \frac{\zeta \left(1+4ϵ\right)}{\zeta \left(1\right)}}`$ (B.53)
$`\times \left(e^\gamma \mu ^2\right)^ϵ\left(4\pi T\right)^{4ϵ}{\displaystyle \frac{T^2}{12}},`$
$`{\displaystyle _𝐩}n_F\left(p\right)p^{12ϵ}`$ $`=`$ $`\left[12^{1+4ϵ}\right]{\displaystyle _𝐩}n_B\left(p\right)p^{12ϵ},`$ (B.54)
This yields

$`2{\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}\mathrm{Re}{\displaystyle _Q}{\displaystyle \frac{1}{Q^2}}{\displaystyle \frac{c^2}{R_0^2+r^2c^2}}_c|_{P_0=ip+\epsilon }`$
$`+{\displaystyle _𝐩}{\displaystyle \frac{n_B\left(p\right)}{p}}c^{1+2ϵ}{\displaystyle _Q}{\displaystyle \frac{1}{Q^2R^2}}|_{P(ip,𝐩/c)}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}{\displaystyle \frac{1}{48}}[{\displaystyle \frac{1}{ϵ^2}}`$
$`(612\mathrm{log}24{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}){\displaystyle \frac{1}{ϵ}}+70.122].`$ (B.55)
For the two terms in (B.51) with a double thermal integral, the averages weighted by $`c^2`$ are given in (C.82) and (C.86). Adding them to (B.55), the final result is

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{1}{P^2Q^2}}{\displaystyle \frac{c^2}{R_0^2+r^2c^2}}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{1}{48}}\right)`$
$`\times \left[{\displaystyle \frac{1}{ϵ^2}}\left(612\mathrm{log}24{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right){\displaystyle \frac{1}{ϵ}}+51.9307\right].`$ (B.56)
Inserting this into (B.51), we obtain the final result (B.47).

The sum-integral (B.48) is evaluated in a similar way to (B.47). Using the representation (A.124) for $`𝒯_R`$, we get

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{q^2}{P^2Q^2r^4}}𝒯_R`$ $`=`$ $`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{q^2}{P^2Q^2r^4}}`$ (B.57)
$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{q^2}{P^2Q^2r^2}}{\displaystyle \frac{c^2}{R_0^2+r^2c^2}}_c.`$
The first sum-integral on the right hand side is given by (B.32). To evaluate the second sum-integral, we apply the sum-integral formula (B.39):

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{q^2}{P^2Q^2r^2\left(R_0^2+r^2c^2\right)}}`$ (B.58)
$`={\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}\mathrm{Re}{\displaystyle _Q}{\displaystyle \frac{p^2+q^2}{Q^2r^2\left(R_0^2+r^2c^2\right)}}|_{P_0=ip+\epsilon }`$
$`+c^{1+2ϵ}{\displaystyle _𝐩}{\displaystyle \frac{n_B\left(p\right)}{p}}p^2{\displaystyle _Q}{\displaystyle \frac{q^2}{Q^2R^2}}|_{P(ip,𝐩/c)}`$
$`+{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_F\left(q\right)}{pq}}{\displaystyle \frac{q^2}{r^2}}\mathrm{Re}{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}`$
$`c^{1+2ϵ}{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_B\left(q\right)}{pq}}{\displaystyle \frac{p^2+r_c^2}{q^2}}`$
$`\times \mathrm{Re}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}.`$
In the terms on the right side with a single thermal integral, the weighted averages over $`c`$ of the integrals over $`Q`$ are given in (C.187), (C.192), and (C.193): After using (B.54) to evaluate the thermal integral, we obtain

$`{\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}\mathrm{Re}{\displaystyle _Q}{\displaystyle \frac{p^2+q^2}{Q^2r^2}}{\displaystyle \frac{c^2}{R_0^2+r^2c^2}}_c|_{P_0=ip+\epsilon }`$
$`+{\displaystyle _𝐩}{\displaystyle \frac{n_B\left(p\right)}{p}}{\displaystyle \frac{1}{p^2}}c^{1+2ϵ}{\displaystyle _Q}{\displaystyle \frac{q^2}{Q^2R^2}}|_{P(ip,𝐩/c)}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{1}{576}}\right)[{\displaystyle \frac{1}{ϵ^2}}`$
$`({\displaystyle \frac{34}{3}}36\mathrm{log}24{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}){\displaystyle \frac{1}{ϵ}}+229.354],`$ (B.59)
For the two terms in (B.58) with a double thermal integral, the averages weighted by $`c^2`$ are given in (C.84), (C.88), and (C.89). Adding them to (B.59), the final result is

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{q^2}{P^2Q^2r^2}}{\displaystyle \frac{c^2}{R_0^2+r^2c^2}}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{1}{576}}\right)`$
$`\times \left[{\displaystyle \frac{1}{ϵ^2}}\left({\displaystyle \frac{118}{3}}52\mathrm{log}24{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right){\displaystyle \frac{1}{ϵ}}+91.002\right].`$
(B.60)
To evaluate (B.49), we use the expression (A.124) for $`𝒯_R`$ and the identity $`PQ=\left(R^2P^2Q^2\right)/2`$ to write it in the form

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{PQ}{P^2Q^2r^4}}𝒯_R={\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{PQ}{P^2Q^2r^4}}`$
$`{\displaystyle _{\left\{P\right\}}}{\displaystyle \frac{1}{P^2}}{\displaystyle _R}{\displaystyle \frac{1}{r^4}}𝒯_R{\displaystyle \frac{1}{2}}c^2_c{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{1}{P^2Q^2r^2}}`$
$`{\displaystyle \frac{1}{2}}{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{1}{P^2Q^2}}{\displaystyle \frac{c^2\left(1c^2\right)}{R_0^2+r^2c^2}}_c.`$ (B.61)
The sum-integrals in the first 3 terms on the right side of (B.61) are given in (B.10), (B.18), (B.31), and (B.34). The last sum-integral before the average weighted by $`c`$ is given in (B.52). The average weighted by $`c^2`$ is given in (B.56). The average weighted by $`c^4`$ can be computed in the same way. In the integrand of the single thermal integral, the weighted averages over $`c`$ of the integrals over $`Q`$ are given in (C.182) and (C.191): After using (B.54) to evaluate the thermal integral, we obtain

$`2{\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}\mathrm{Re}{\displaystyle _Q}{\displaystyle \frac{1}{Q^2}}{\displaystyle \frac{c^4}{R_0^2+r^2c^2}}_c|_{P_0=ip+\epsilon }`$
$`+{\displaystyle _𝐩}{\displaystyle \frac{n_B\left(p\right)}{p}}c^{1+2ϵ}{\displaystyle _Q}{\displaystyle \frac{1}{Q^2R^2}}|_{P(ip,𝐩/c)}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left[\left({\displaystyle \frac{7}{72}}{\displaystyle \frac{1}{6}}\mathrm{log}2\right){\displaystyle \frac{1}{ϵ}}+0.2150\right],`$
(B.62)
For the two terms with a double thermal integral, the averages weighted by $`c^4`$ are given in (C.83) and (C.87). Adding them to (B.62), we obtain

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{1}{P^2Q^2}}{\displaystyle \frac{c^4}{R_0^2+r^2c^2}}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$
$`\times \left[\left({\displaystyle \frac{7}{72}}{\displaystyle \frac{1}{6}}\mathrm{log}2\right){\displaystyle \frac{1}{ϵ}}+0.1359\right],`$ (B.63)
We finally need to evaluate (B.50). Applying (B.39) gives

$`{\displaystyle _{\left\{PQ\right\}}}{\displaystyle \frac{r^2p^2}{P^2q^2Q_0^2R^2}}𝒯_Q`$ $`=`$ $`\left[{\displaystyle _𝐩}{\displaystyle \frac{n_B\left(p\right)}{p}}+{\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}\right]`$ (B.64)
$`\times \mathrm{Re}{\displaystyle _Q}{\displaystyle \frac{p^2q^2}{Q^2r^2\left(R_0^2+r^2c^2\right)}}_c|_{P_0=ip}`$
$`+{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_F\left(q\right)}{pq}}\mathrm{Re}{\displaystyle \frac{r_c^2p^2}{q^2}}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+iϵ,q,r_c)}}c^{1+2ϵ}_c`$
$`+{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_B\left(p\right)n_F\left(q\right)}{pq}}\mathrm{Re}{\displaystyle \frac{r_c^2p^2}{q^2}}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+iϵ,q,r_c)}}c^{1+2ϵ}_c`$
$`+{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_B\left(q\right)}{pq}}\mathrm{Re}{\displaystyle \frac{p^2q^2}{r^2}}{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+iϵ,q,rc)}}_c.`$
In the terms on the right side, with a single thermal factor, the weighted average is given in Eq. (C.194), After using Eq. (B.54) to evaluate the thermal integral, we obtain

$`\left[{\displaystyle _𝐩}{\displaystyle \frac{n_B\left(p\right)}{p}}+{\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}\right]{\displaystyle _Q}{\displaystyle \frac{p^2q^2}{Q^2r^2\left(R_0^2+r^2c^2\right)}}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\pi ^2}{24}}\right).`$ (B.65)
The terms with two thermal factors are given in Eqs. (C.85), (C.90) and (C.91). Adding them to (B.65), we finally obtain (B.50).

## Appendix C Integrals

Dimensional regularization can be used to regularize both the ultraviolet divergences and infrared divergences in 3-dimensional integrals over momenta. The spacial dimension is generalized to $`d=32ϵ`$ 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

$$_𝐩\left(\frac{e^\gamma \mu ^2}{4\pi }\right)^ϵ\frac{d^{32ϵ}p}{\left(2\pi \right)^{32ϵ}}.$$
(C.66)
### C.1 3-dimensional integrals

We require one integral that does not involve the Bose-Einstein distribution function. The momentum scale in these integrals is set by the mass $`m=m_D`$. The one-loop integral is

$`{\displaystyle _𝐩}{\displaystyle \frac{1}{p^2+m^2}}`$ $`=`$ $`{\displaystyle \frac{m}{4\pi }}\left({\displaystyle \frac{\mu }{2m}}\right)^{2ϵ}\left[1+2ϵ\right].`$ (C.67)
The error is one order higher in $`ϵ`$ than the smallest term shown.

### C.2 Thermal integrals

The thermal integrals involve the Fermi-Dirac distribution $`n_F\left(p\right)=1/\left(e^{\beta p}+1\right)`$. The one-loop integrals can all be obtained from the general formula

$`{\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}p^{2\alpha }`$ $`=`$ $`\left(12^{12\alpha +2ϵ}\right){\displaystyle \frac{\zeta \left(2+2\alpha 2ϵ\right)}{4\pi ^2}}`$ (C.68)
$`\times {\displaystyle \frac{\mathrm{\Gamma }\left(2+2\alpha 2ϵ\right)\mathrm{\Gamma }\left(\frac{1}{2}\right)}{\mathrm{\Gamma }\left(\frac{3}{2}ϵ\right)}}\left(e^\gamma \mu ^2\right)^ϵT^{2+2\alpha 2ϵ}.`$
The simple two-loop thermal integrals are

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_F\left(q\right)}{pq}}{\displaystyle \frac{1}{r^2}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}{\displaystyle \frac{1}{3}}\left[1\mathrm{log}2\right],`$ (C.69)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_F\left(q\right)}{pq}}{\displaystyle \frac{q^2}{r^4}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$ (C.70)
$`\times \left({\displaystyle \frac{1}{36}}\right)\left[5+6\gamma +6\mathrm{log}26{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}+3.076ϵ\right],`$
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_B\left(p\right)n_F\left(q\right)}{pq}}{\displaystyle \frac{p^2}{r^4}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$ (C.71)
$`\times \left({\displaystyle \frac{1}{36}}\right)\left[76\gamma 18\mathrm{log}2+6{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}+29.508ϵ\right],`$
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_B\left(p\right)n_F\left(q\right)}{pq}}{\displaystyle \frac{q^2}{r^4}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$ (C.72)
$`\times \left({\displaystyle \frac{1}{18}}\right)\left[16\gamma 12\mathrm{log}2+6{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}+31.134ϵ\right].`$
We also need some more complicated 2-loop thermal integrals that involve the triangle function defined in Eq. (B.46):

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{r^4}{q^2\mathrm{\Delta }(p,q,r)}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{7}{96}}\right)[{\displaystyle \frac{1}{ϵ^2}}+\left({\displaystyle \frac{22}{7}}+2\gamma +2\mathrm{log}2+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}{\displaystyle \frac{7}{20}}\zeta \left(3\right)\right){\displaystyle \frac{1}{ϵ}}+47.2406],`$ (C.73)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{r^2}{\mathrm{\Delta }(p,q,r)}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}({\displaystyle \frac{1}{48}})[{\displaystyle \frac{1}{ϵ^2}}+2(1+\gamma +\mathrm{log}2+{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}){\displaystyle \frac{1}{ϵ}}`$ (C.74)
$`+4+4\gamma +{\displaystyle \frac{\pi ^2}{2}}+4\gamma \mathrm{log}26\mathrm{log}^22+4\mathrm{log}24\gamma _1+4(1+\gamma +\mathrm{log}2){\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}+2{\displaystyle \frac{\zeta ^{\prime \prime }\left(1\right)}{\zeta \left(1\right)}}],`$
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{p^4}{q^2\mathrm{\Delta }(p,q,r)}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{49\zeta \left(3\right)}{1920}}\right)\left[{\displaystyle \frac{1}{ϵ}}+2+2\mathrm{log}2+2{\displaystyle \frac{\zeta ^{}\left(3\right)}{\zeta \left(3\right)}}+2{\displaystyle \frac{\zeta ^{}\left(3\right)}{\zeta \left(3\right)}}\right],`$ (C.75)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{p^2\left(p^2+q^2\right)}{r^2\mathrm{\Delta }(p,q,r)}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{1}{96}}\right)\left[{\displaystyle \frac{1}{ϵ^2}}+\left({\displaystyle \frac{26}{3}}+10\gamma 6{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}+10\mathrm{log}2\right){\displaystyle \frac{1}{ϵ}}+41.1586\right],`$ (C.76)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{p^2}{\mathrm{\Delta }(p,q,r)}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{1}{96}}\right)\left[{\displaystyle \frac{1}{ϵ^2}}+2\left(1+\gamma +\mathrm{log}2+{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right){\displaystyle \frac{1}{ϵ}}+37.0573\right],`$ (C.77)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_B\left(q\right)}{q}}{\displaystyle \frac{p^2}{\mathrm{\Delta }(p,q,r)}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{1}{96}}\right)\left[{\displaystyle \frac{1}{ϵ^2}}+2\left(1+\gamma \mathrm{log}2+{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right){\displaystyle \frac{1}{ϵ}}+19.2257\right],`$ (C.78)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_B\left(q\right)}{q}}{\displaystyle \frac{p^4}{q^2\mathrm{\Delta }(p,q,r)}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{7\zeta \left(3\right)}{1920}}\right)\left[{\displaystyle \frac{1}{ϵ}}+2{\displaystyle \frac{2}{7}}\mathrm{log}2+2{\displaystyle \frac{\zeta ^{}\left(3\right)}{\zeta \left(3\right)}}+2{\displaystyle \frac{\zeta ^{}\left(3\right)}{\zeta \left(3\right)}}\right],`$ (C.79)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_B\left(q\right)}{q}}{\displaystyle \frac{r^4}{q^2\mathrm{\Delta }(p,q,r)}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{1}{24}}\right)\left[{\displaystyle \frac{1}{ϵ^2}}+\left(4+2\gamma 5\mathrm{log}2{\displaystyle \frac{7\zeta \left(3\right)}{80}}+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right){\displaystyle \frac{1}{ϵ}}+18.1551\right],`$ (C.80)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_B\left(q\right)}{q}}{\displaystyle \frac{r^2}{\mathrm{\Delta }(p,q,r)}}`$ $`=`$ $`{\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{1}{96}}\right)\left[{\displaystyle \frac{1}{ϵ^2}}+2\left(1+\gamma +5\mathrm{log}2+{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right){\displaystyle \frac{1}{ϵ}}+84.2513\right].`$ (C.81)
The most difficult thermal integrals to evaluate involve both the triangle function and the HTL average defined in (A.125). There are 2 sets of these integrals. The first set is

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_F\left(q\right)}{pq}}\text{Re}c^2{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left[1.458\times 10^2\right],`$ (C.82)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_F\left(q\right)}{pq}}\text{Re}c^4{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left[1.7715\times 10^2\right],`$ (C.83)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_F\left(q\right)}{pq}}\text{Re}{\displaystyle \frac{q^2}{r^2}}c^2{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left[1.1578\times 10^2\right],`$ (C.84)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_B\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}\text{Re}{\displaystyle \frac{p^2q^2}{r^2}}{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+iϵ,q,rc)}}`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left[0.17811\right].`$ (C.85)
The second set of these integrals involve the variable $`r_c=\left|𝐩+𝐪/c\right|`$:

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_B\left(q\right)}{pq}}\mathrm{Re}c^{1+2ϵ}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left[0.19678\right],`$ (C.86)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_B\left(q\right)}{pq}}\mathrm{Re}c^{1+2ϵ}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left[4.8368\times 10^2\right],`$ (C.87)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_B\left(q\right)}{pq}}{\displaystyle \frac{p^2}{q^2}}\mathrm{Re}c^{1+2ϵ}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}{\displaystyle \frac{1}{96}}\left[{\displaystyle \frac{1}{ϵ}}+7.7702\right],`$ (C.88)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_B\left(q\right)}{pq}}\mathrm{Re}c^{1+2ϵ}{\displaystyle \frac{r_c^2}{q^2}}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}{\displaystyle \frac{118\mathrm{log}2}{288}}\left[{\displaystyle \frac{1}{ϵ}}+7.79693\right],`$ (C.89)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(p\right)}{q}}\text{Re}{\displaystyle \frac{r_c^2p^2}{q^2}}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+iϵ,q,r_c)}}c^{1+2ϵ}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{1}{24}}\right)[{\displaystyle \frac{1}{ϵ^2}}+(2+2\gamma +2\mathrm{log}2+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}){\displaystyle \frac{1}{ϵ}}`$
$`+40.316],`$ (C.90)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_B\left(p\right)}{p}}{\displaystyle \frac{n_F\left(p\right)}{q}}\text{Re}{\displaystyle \frac{r_c^2p^2}{q^2}}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+iϵ,q,r_c)}}c^{1+2ϵ}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{1}{12}}\right)[{\displaystyle \frac{1}{ϵ^2}}+(2+2\gamma +4\mathrm{log}2+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}){\displaystyle \frac{1}{ϵ}}`$
$`+52.953].`$ (C.91)
The simplest way to evaluate integrals like (C.69)–(C.72) whose integrands factor into separate functions of $`p`$, $`q`$, and $`r`$ is to Fourier transform to coordinate space where they reduce to an integral over a single coordinate $`𝐑`$:

$`{\displaystyle _{\mathrm{𝐩𝐪}}}f\left(p\right)g\left(q\right)h\left(r\right)`$ $`=`$ $`{\displaystyle _𝐑}\stackrel{~}{f}\left(R\right)\stackrel{~}{g}\left(R\right)\stackrel{~}{h}\left(R\right).`$ (C.92)
The Fourier transform is

$`\stackrel{~}{f}\left(R\right)={\displaystyle _𝐩}e^{i𝐩𝐑}f\left(p\right),`$ (C.93)
and the dimensionally regularized coordinate integral is

$`{\displaystyle _𝐑}`$ $`=`$ $`\left({\displaystyle \frac{e^\gamma \mu ^2}{4\pi }}\right)^ϵ{\displaystyle d^{32ϵ}R}.`$ (C.94)
The Fourier transforms we need are

$`{\displaystyle _𝐩}p^{2\alpha }e^{i𝐩𝐑}`$ $`=`$ $`{\displaystyle \frac{1}{8\pi }}{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3}{2}+\alpha ϵ\right)}{\mathrm{\Gamma }\left(\frac{1}{2}\right)\mathrm{\Gamma }\left(\alpha \right)}}\left(e^\gamma \mu ^2\right)^ϵ`$ (C.95)
$`\times \left({\displaystyle \frac{2}{R}}\right)^{3+2\alpha 2ϵ},`$
$`{\displaystyle _𝐩}{\displaystyle \frac{n\left(p\right)}{p}}p^{2\alpha }e^{i𝐩𝐑}`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle \frac{1}{\mathrm{\Gamma }\left(\frac{1}{2}\right)}}\left(e^\gamma \mu ^2\right)^ϵ\left({\displaystyle \frac{2}{R}}\right)^{\frac{1}{2}ϵ}`$ (C.96)
$`\times {\displaystyle _0^{\mathrm{}}}dpp^{2\alpha +\frac{1}{2}ϵ}n\left(p\right)J_{\frac{1}{2}ϵ}\left(pR\right).`$
If $`\alpha `$ is an even integer, the Fourier transform (C.96) is particularly simple in the limit $`d3`$:

$`{\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}e^{i𝐩𝐑}={\displaystyle \frac{T}{4\pi R}}\left({\displaystyle \frac{1}{x}}\text{csch}x\right),`$ (C.97)
$`{\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}p^2e^{i𝐩𝐑}={\displaystyle \frac{\pi T^3}{2R}}\left(\text{csch}^3x+{\displaystyle \frac{1}{2}}\text{csch}x{\displaystyle \frac{1}{x^3}}\right),`$ (C.98)
where $`x=\pi RT`$

We can use these simple expressions only if the integral over the coordinate $`𝐑`$ in (C.92) converges for $`d=3`$. Otherwise, we must first make subtractions on the integrand to make the integral convergent.

The integrals (C.69)–(C.72) can be evaluated directly by applying the Fourier transform formula (C.92) in the limit $`ϵ0`$.

The integrals (C.73)–(C.75) can be evaluated by first averaging over angles. The triangle function can be expressed as

$$\mathrm{\Delta }(p,q,r)=4p^2q^2\left(1\mathrm{cos}^2\theta \right),$$
(C.99)
where $`\theta `$ is the angle between $`𝐩`$ and $`𝐪`$. For example, the angle average for (C.73) is

$`{\displaystyle \frac{r^4}{\mathrm{\Delta }(p,q,r)}}_{\widehat{𝐩}\widehat{𝐪}}`$ $`=`$ $`{\displaystyle \frac{w\left(ϵ\right)}{8}}{\displaystyle _1^{+1}}𝑑x\left(1x^2\right)^{1ϵ}`$ (C.100)
$`\times {\displaystyle \frac{\left(p^2+q^2+2pqx\right)^2}{p^2q^2}}.`$
After integrating over $`x`$ and inserting the result into (C.73), the integral reduces to

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{r^4}{q^2\mathrm{\Delta }(p,q,r)}}`$
$`={\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}\left({\displaystyle \frac{12ϵ}{8ϵ}}{\displaystyle \frac{p^2}{q^4}}+{\displaystyle \frac{76ϵ}{8ϵ}}{\displaystyle \frac{1}{q^2}}\right).`$ (C.101)
The integrals over $`𝐩`$ and $`𝐪`$ factor into separate integrals that can be evaluated using (C.68). After averaging over angles, the integrals (C.74) and (C.75) reduce to

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{r^2}{\mathrm{\Delta }(p,q,r)}}`$
$`={\displaystyle \frac{12ϵ}{4ϵ}}{\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle _𝐪}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{1}{q^2}},`$ (C.102)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{p^4}{q^2\mathrm{\Delta }(p,q,r)}}`$
$`={\displaystyle \frac{12ϵ}{8ϵ}}{\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}p^2{\displaystyle _𝐪}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{1}{q^4}}.`$ (C.103)
The integral (C.76) can be evaluated by using the identity

$$\frac{p^2+q^2}{r^2\mathrm{\Delta }(p,q,r)}_{\widehat{𝐩}\widehat{𝐪}}=\frac{1}{2ϵ}\frac{1}{r^4}_{\widehat{𝐩}\widehat{𝐪}}+\frac{12ϵ}{8ϵ}\frac{1}{p^2q^2}.$$
(C.104)
The identity can be proved by expressing the angular averages in terms of integrals over the cosine of the angle between $`𝐩`$ and $`𝐪`$ as in (C.100), and then integrating by parts. Inserting the identity (C.104) into (C.76), the integral reduces to

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{p^2\left(p^2+q^2\right)}{r^2\mathrm{\Delta }(p,q,r)}}`$ $`=`$ $`{\displaystyle \frac{1}{2ϵ}}{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{p^2}{r^4}}`$ (C.105)
$`+{\displaystyle \frac{12ϵ}{8ϵ}}{\displaystyle _𝐩}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle _𝐪}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{1}{q^2}}.`$
The integral in the first term on the right is given in (C.70), while the second term can be evaluated using (C.68).

The integral (C.82) can be evaluated directly in three dimensions by first averaging over $`c`$ and $`x`$, and then integrate the resulting functions numerically over $`p`$ and $`q`$.

To evaluate the weighted averages over $`c`$ of the thermal integrals in Eqs.(C.83)– (C.85), we first isolate the divergent parts, which come from the region $`pq0`$. We write the product of thermal functions in the form

$`n_F\left(p\right)n_F\left(q\right)`$ $`=`$ $`\left(n_F\left(p\right)n_F\left(q\right){\displaystyle \frac{s^2n_F^2\left(s\right)}{pq}}\right)`$ (C.106)
$`+{\displaystyle \frac{s^2n_F^2\left(s\right)}{pq}},`$
where $`s=\left(p+q\right)/2`$. In the difference term, the HTL average over $`c`$ and the angular average over $`x=\widehat{𝐩}\widehat{𝐪}`$ can be calculated in three dimensions:

$`\mathrm{Re}c^4{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}_{c,x}`$ $`=`$ $`{\displaystyle \frac{2\left(p^2+q^2\right)}{3\left(p^2q^2\right)^2}}`$
$`+{\displaystyle \frac{1}{12pq}}\mathrm{log}{\displaystyle \frac{p+q}{\left|pq\right|}}{\displaystyle \frac{\left(3p^2+q^2\right)\left(p^2+3q^2\right)}{6\left(p^2q^2\right)^3}}\mathrm{log}\left(p/q\right),`$
$`\mathrm{Re}c^2{\displaystyle \frac{q^2}{r^2}}{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}_{c,x}`$ $`=`$ $`{\displaystyle \frac{q^2}{3\left(p^2q^2\right)^2}}`$ (C.108)
$`\times [2{\displaystyle \frac{1}{2}}\mathrm{log}{\displaystyle \frac{\left|p^2q^2\right|}{pq}}`$
$`{\displaystyle \frac{p^2+q^2}{4pq}}\mathrm{log}{\displaystyle \frac{p+q}{\left|pq\right|}}{\displaystyle \frac{p^2+q^2}{p^2q^2}}\mathrm{log}(p/q)],`$
$`\text{Re}{\displaystyle \frac{p^2q^2}{r^2}}{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+iϵ,q,rc)}}_{c,x}`$ $`=`$ $`{\displaystyle \frac{1}{4pq\left(p^2q^2\right)}}`$ (C.109)
$`\times \left[\left(p^2+q^2\right)\mathrm{log}{\displaystyle \frac{p+q}{\left|pq\right|}}2pq\mathrm{log}{\displaystyle \frac{\left|p^2q^2\right|}{pq}}\right].`$
The remaining 2-dimensional integral over $`p`$ and $`q`$ can then be evaluated numerically:

$`{\displaystyle _{\mathrm{𝐩𝐪}}}\left({\displaystyle \frac{n_F\left(p\right)n_F\left(q\right)}{pq}}{\displaystyle \frac{s^2n_F^2\left(s\right)}{p^2q^2}}\right)`$
$`\times \mathrm{Re}c^4{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}[8.980\times 10^3],`$ (C.110)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}\left({\displaystyle \frac{n_F\left(p\right)n_F\left(q\right)}{pq}}{\displaystyle \frac{s^2n_F^2\left(s\right)}{p^2q^2}}\right){\displaystyle \frac{q^2}{r^2}}`$
$`\times \mathrm{Re}c^2{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}[7.792\times 10^3],`$ (C.111)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}\left({\displaystyle \frac{n_B\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{s^2n_B\left(s\right)n_F\left(s\right)}{p^2q^2}}\right)`$
$`\times \text{Re}{\displaystyle \frac{p^2q^2}{r^2}}{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+iϵ,q,rc)}}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left[0.17811\right].`$ (C.112)
The integrals involving the $`n_F^2\left(s\right)`$ term in (C.106) are divergent, so the HTL average over $`c`$ and the angular average over $`x=\widehat{𝐩}\widehat{𝐪}`$ must be calculated in $`32ϵ`$ dimensions. The first step in the calculation of the $`n^2\left(s\right)`$ term is to change variables from $`𝐩`$ and $`𝐪`$ to $`s=\left(p+q\right)/2`$, $`\beta =4pq/\left(p+q\right)^2`$, and $`x=\widehat{𝐩}\widehat{𝐪}`$:

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{s^2n_F^2\left(s\right)}{p^2q^2}}f(p,q,r)={\displaystyle \frac{64}{\left(4\pi \right)^4}}\left[\left(e^\gamma \mu ^2\right)^ϵ{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{3}{2}\right)}{\mathrm{\Gamma }\left(\frac{3}{2}ϵ\right)}}\right]^2`$
$`\times {\displaystyle _0^{\mathrm{}}}dss^{14ϵ}n_F^2\left(s\right)s^2{\displaystyle _0^1}d\beta \beta ^{2ϵ}(1\beta )^{1/2}`$
$`\times f(s_+,s_{},r)+f(s_{},s_+,r)_x,`$ (C.113)
where $`s_\pm =s\left[1\pm \sqrt{1\beta }\right]`$ and $`r=s\left[42\beta \left(1x\right)\right]^{1/2}`$. The 2 terms inside the average over $`x`$ come from the regions $`p>q`$ and $`p<q`$, respectively. The integral over $`s`$ is easily evaluated:

$`{\displaystyle _0^{\mathrm{}}}dss^{14ϵ}n_F^2\left(s\right)=\mathrm{\Gamma }(24ϵ)[(12^{4ϵ})\zeta (14ϵ)`$
$`+(12^{1+4ϵ})\zeta (24ϵ)]T^{24ϵ},`$ (C.114)
$`{\displaystyle _0^{\mathrm{}}}𝑑ss^{14ϵ}n_F\left(s\right)n_B\left(s\right)`$
$`=2^{2+4ϵ}\mathrm{\Gamma }\left(24ϵ\right)\zeta \left(24ϵ\right)T^{24ϵ}.`$ (C.115)
It remains only to evaluate the averages over $`c`$ and $`x`$ and the integral over $`\beta `$.

The first step in the calculation of the $`n_F^2\left(s\right)`$ term of (C.83) is to decompose the integrand into 2 terms:

$$\frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}=\frac{1}{2}\underset{\pm }{}\frac{1}{\left(p+i\epsilon \pm q\right)^2r^2c^2}.$$
(C.116)
The weighted averages over $`c`$ gives a hypergeometric function:

$`{\displaystyle \frac{c^4}{\left(p+i\epsilon \pm q\right)^2r^2c^2}}_c={\displaystyle \frac{3}{\left(32ϵ\right)\left(52ϵ\right)}}`$
$`\times {\displaystyle \frac{1}{\left(p+i\epsilon \pm q\right)^2}}F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{5}{2},1}{\frac{7}{2}ϵ}}|{\displaystyle \frac{r^2}{\left(p+i\epsilon \pm q\right)^2}}\right).`$ (C.117)
In the $`+q`$ case of (C.117), the $`i\epsilon `$ prescription is unnecessary. The argument of the hypergeometric function can be written $`1\beta y`$, where $`y=\left(1x\right)/2`$. After using a transformation formula to change the argument to $`\beta y`$, we can evaluate the angular average over $`x`$ to obtain hypergeometric functions with argument $`\beta `$. For example, the average over $`x`$ of (C.117) is

$`F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{5}{2},1}{\frac{7}{2}ϵ}}\right|{\displaystyle \frac{r^2}{\left(p+q\right)^2}})_x={\displaystyle \frac{52ϵ}{2ϵ}}[F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,\frac{5}{2},1}{22ϵ,1+ϵ}}\right|\beta )`$
$`{\displaystyle \frac{\left(1\right)_ϵ\left(1\right)_{2ϵ}\left(2\right)_{2ϵ}\left(\frac{5}{2}\right)_ϵ}{\left(1\right)_ϵ\left(2\right)_{3ϵ}}}\beta ^ϵF\left({\displaystyle \genfrac{}{}{0pt}{}{12ϵ,\frac{5}{2}ϵ}{23ϵ}}\right|\beta )],`$ (C.118)
where $`\left(a\right)_b`$ is Pochhammer’s symbol which is defined in (C.212). Integrating over $`\beta `$, we obtain hypergeometric functions with argument 1:

$`s^2{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}{\displaystyle \frac{c^4}{\left(p+q\right)^2r^2c^2}}_{c,x}={\displaystyle \frac{1}{4ϵ}}`$
$`\times {\displaystyle \frac{\left(1\right)_ϵ\left(2\right)_{2ϵ}}{\left(1\right)_ϵ}}[{\displaystyle \frac{\left(1\right)_{2ϵ}\left(1\right)_ϵ}{\left(\frac{5}{2}\right)_{2ϵ}\left(2\right)_{2ϵ}\left(1\right)_ϵ}}F\left({\displaystyle \genfrac{}{}{0pt}{}{12ϵ,1ϵ,\frac{5}{2},1}{\frac{5}{2}2ϵ,22ϵ,1+ϵ}}\right|1)`$
$`{\displaystyle \frac{\left(1\right)_{3ϵ}\left(1\right)_{2ϵ}\left(\frac{5}{2}\right)_ϵ}{\left(\frac{5}{2}\right)_{3ϵ}\left(2\right)_{3ϵ}}}F\left({\displaystyle \genfrac{}{}{0pt}{}{13ϵ,12ϵ,\frac{5}{2}ϵ}{\frac{5}{2}3ϵ,23ϵ}}\right|1)].`$ (C.119)
Expanding in powers of $`ϵ`$, we obtain

$`s^2{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}{\displaystyle \frac{c^4}{\left(p+q\right)^2r^2c^2}}_{c,x}`$
$`={\displaystyle \frac{\pi ^2}{72}}\left(1+10.8408ϵ\right).`$ (C.120)
In the $`q`$ case of (C.117), the argument of the hypergeometric functions can be written $`\left(1\beta y\right)/\left(1\beta \pm i\epsilon \right)`$, where $`y=\left(1x\right)/2`$ and the prescriptions $`+i\epsilon `$ and $`i\epsilon `$ correspond to the regions $`p>q`$ and $`p<q`$, respectively. These regions correspond to the two terms inside the average over $`x`$ in (C.113). In order to obtain an analytic result in terms of hypergeometric functions, it is necessary to integrate over $`\beta `$ before averaging over $`x`$. The integrals over $`\beta `$ can be evaluated by first using a transformation formula to change the argument of the hypergeometric function to $`\beta \left(1y\right)/\left(1\beta \right)`$ and then using the integration formula (C.219) to obtain hypergeometric functions with arguments $`y`$ or $`1y`$:

$`{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{3/2}F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{3}{2},1}{\frac{5}{2}ϵ}}|{\displaystyle \frac{1\beta y}{1\beta +i\epsilon }}\right)`$
$`={\displaystyle \frac{32ϵ}{ϵ}}{\displaystyle \frac{\left(1\right)_{2ϵ}}{\left(\frac{1}{2}\right)_{2ϵ}}}F\left({\displaystyle \genfrac{}{}{0pt}{}{12ϵ,1}{1+ϵ}}|1y\right)`$
$`{\displaystyle \frac{32ϵ}{ϵ}}{\displaystyle \frac{\left(1\right)_ϵ}{\left(\frac{1}{2}\right)_ϵ}}\left(1y\right)^{1/2}F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{1}{2}2ϵ,1}{\frac{1}{2}+ϵ}}|1y\right)`$
$`+{\displaystyle \frac{3}{2ϵ\left(13ϵ\right)}}e^{i\pi ϵ}\left(1\right)_ϵ\left(\frac{5}{2}\right)_ϵ\left(1y\right)^ϵF\left({\displaystyle \genfrac{}{}{0pt}{}{13ϵ,\frac{3}{2}ϵ}{23ϵ}}|y\right).`$
After averaging over $`x`$, we obtain hypergeometric functions with argument 1:

$`s^2{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}{\displaystyle \frac{c^2}{\left(p+i\epsilon q\right)^2r^2c^2}}_{c,x}`$
$`={\displaystyle \frac{1}{4ϵ}}{\displaystyle \frac{\left(1\right)_{2ϵ}}{\left(\frac{1}{2}\right)_{2ϵ}}}F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,12ϵ,1}{22ϵ,1+ϵ}}|1\right)`$
$`{\displaystyle \frac{1}{2ϵ}}{\displaystyle \frac{\left(2\right)_{2ϵ}\left(1\right)_ϵ\left(\frac{1}{2}\right)_ϵ}{\left(1\right)_ϵ\left(\frac{1}{2}\right)_ϵ\left(\frac{3}{2}\right)_{2ϵ}}}F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{1}{2}ϵ,\frac{1}{2}2ϵ,1}{\frac{3}{2}2ϵ,\frac{1}{2}+ϵ}}|1\right)`$
$`+{\displaystyle \frac{1}{8ϵ\left(13ϵ\right)}}e^{i\pi ϵ}{\displaystyle \frac{\left(2\right)_{2ϵ}\left(1\right)_{2ϵ}\left(1\right)_ϵ\left(\frac{3}{2}\right)_ϵ}{\left(1\right)_ϵ\left(2\right)_{3ϵ}}}`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,13ϵ,\frac{3}{2}ϵ}{23ϵ,23ϵ}}|1\right).`$ (C.122)
Expanding in powers of $`ϵ`$ and then taking the real parts, we obtain

$`\mathrm{Re}s^2{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}{\displaystyle \frac{c^4}{\left(p+i\epsilon q\right)^2r^2c^2}}_{c,x}`$
$`={\displaystyle \frac{12+\pi ^2}{72}}\left(1+1.10518ϵ\right).`$ (C.123)
To evaluate the subtraction in the integral (C.111), we use the identity $`q^2=\left(r^2+q^2p^22𝐩𝐪\right)/2`$. The integral with $`q^2p^2`$ in the numerator is purely imaginary. Thus the real part of the integral can be expressed as

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{s^2n^2\left(s\right)}{p^2q^2}}{\displaystyle \frac{q^2}{r^2}}\mathrm{Re}c^2{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}_c`$
$`={\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{s^2n^2\left(s\right)}{p^2q^2}}\left({\displaystyle \frac{1}{2}}{\displaystyle \frac{𝐩𝐪}{r^2}}\right)\mathrm{Re}c^2{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}_c.`$
(C.124)
The first term in Eq. (C.124) is decomposed into 2 terms:

$$\frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}=\frac{1}{2}\underset{\pm }{}\frac{1}{\left(p+i\epsilon \pm q\right)^2r^2c^2}.$$
(C.125)
The weighted averages over $`c`$ give hypergeometric functions:

$`{\displaystyle \frac{c^2}{\left(p+i\epsilon \pm q\right)^2r^2c^2}}_c`$
$`={\displaystyle \frac{1}{32ϵ}}{\displaystyle \frac{1}{\left(p+i\epsilon \pm q\right)^2}}F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{3}{2},1}{\frac{5}{2}ϵ}}|{\displaystyle \frac{r^2}{\left(p+i\epsilon \pm q\right)^2}}\right),`$ (C.126)
$`{\displaystyle \frac{c^4}{\left(p+i\epsilon \pm q\right)^2r^2c^2}}_c={\displaystyle \frac{3}{\left(32ϵ\right)\left(52ϵ\right)}}`$
$`\times {\displaystyle \frac{1}{\left(p+i\epsilon \pm q\right)^2}}F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{5}{2},1}{\frac{7}{2}ϵ}}|{\displaystyle \frac{r^2}{\left(p+i\epsilon \pm q\right)^2}}\right).`$ (C.127)
In the $`+q`$ case of (C.126), the $`i\epsilon `$ prescription is unnecessary. The argument of the hypergeometric function can be written $`1\beta y`$, where $`y=\left(1x\right)/2`$. After using a transformation formula to change the argument to $`\beta y`$, we can evaluate the angular average over $`x`$ to obtain hypergeometric functions with argument $`\beta `$. For example, the average over $`x`$ of (C.126) is

$`F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{3}{2},1}{\frac{5}{2}ϵ}}\right|{\displaystyle \frac{r^2}{\left(p+q\right)^2}})_x={\displaystyle \frac{32ϵ}{2ϵ}}[F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,\frac{3}{2},1}{22ϵ,1+ϵ}}\right|\beta )`$
$`{\displaystyle \frac{\left(1\right)_ϵ\left(1\right)_{2ϵ}\left(2\right)_{2ϵ}\left(\frac{3}{2}\right)_ϵ}{\left(1\right)_ϵ\left(2\right)_{3ϵ}}}\beta ^ϵF\left({\displaystyle \genfrac{}{}{0pt}{}{12ϵ,\frac{3}{2}ϵ}{23ϵ}}\right|\beta )],`$ (C.128)
where $`\left(a\right)_b`$ is Pochhammer’s symbol which is defined in (C.212). Integrating over $`\beta `$, we obtain hypergeometric functions with argument 1:

$`s^2{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}{\displaystyle \frac{c^2}{\left(p+q\right)^2r^2c^2}}_{c,x}`$ (C.129)
$`={\displaystyle \frac{1}{4ϵ}}{\displaystyle \frac{\left(1\right)_ϵ\left(2\right)_{2ϵ}}{\left(1\right)_ϵ}}[{\displaystyle \frac{\left(1\right)_{2ϵ}\left(1\right)_ϵ}{\left(\frac{3}{2}\right)_{2ϵ}\left(2\right)_{2ϵ}\left(1\right)_ϵ}}`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{12ϵ,1ϵ,\frac{3}{2},1}{\frac{3}{2}2ϵ,22ϵ,1+ϵ}}|1\right){\displaystyle \frac{\left(1\right)_{3ϵ}\left(1\right)_{2ϵ}\left(\frac{3}{2}\right)_ϵ}{\left(\frac{3}{2}\right)_{3ϵ}\left(2\right)_{3ϵ}}}`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{13ϵ,12ϵ,\frac{3}{2}ϵ}{\frac{3}{2}3ϵ,23ϵ}}\right|1)].`$ (C.130)
Expanding in powers of $`ϵ`$, we obtain

$`s^2{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}{\displaystyle \frac{c^2}{\left(p+q\right)^2r^2c^2}}_{c,x}`$
$`={\displaystyle \frac{\pi ^2}{24}}+O\left(ϵ\right).`$ (C.131)
In the $`q`$ case of (C.126), the argument of the hypergeometric functions can be written $`\left(1\beta y\right)/\left(1\beta \pm i\epsilon \right)`$, where $`y=\left(1x\right)/2`$ and the prescriptions $`+i\epsilon `$ and $`i\epsilon `$ correspond to the regions $`p>q`$ and $`p<q`$, respectively. These regions correspond to the two terms inside the average over $`x`$ in (C.113). In order to obtain an analytic result in terms of hypergeometric functions, it is necessary to integrate over $`\beta `$ before averaging over $`x`$. The integrals over $`\beta `$ can be evaluated by first using a transformation formula to change the argument of the hypergeometric function to $`\beta \left(1y\right)/\left(1\beta \right)`$ and then using the integration formula (C.219) to obtain hypergeometric functions with arguments $`y`$ or $`1y`$:

$`{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{3/2}F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{3}{2},1}{\frac{5}{2}ϵ}}|{\displaystyle \frac{1\beta y}{1\beta +i\epsilon }}\right)`$
$`={\displaystyle \frac{32ϵ}{ϵ}}{\displaystyle \frac{\left(1\right)_{2ϵ}}{\left(\frac{1}{2}\right)_{2ϵ}}}F\left({\displaystyle \genfrac{}{}{0pt}{}{12ϵ,1}{1+ϵ}}|1y\right)`$
$`{\displaystyle \frac{32ϵ}{ϵ}}{\displaystyle \frac{\left(1\right)_ϵ}{\left(\frac{1}{2}\right)_ϵ}}\left(1y\right)^{1/2}F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{1}{2}2ϵ,1}{\frac{1}{2}+ϵ}}|1y\right)`$
$`+{\displaystyle \frac{3}{2ϵ\left(13ϵ\right)}}e^{i\pi ϵ}\left(1\right)_ϵ\left(\frac{5}{2}\right)_ϵ\left(1y\right)^ϵ`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{13ϵ,\frac{3}{2}ϵ}{23ϵ}}|y\right).`$ (C.132)
After averaging over $`x`$, we obtain hypergeometric functions with argument 1:

$`s^2{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}{\displaystyle \frac{c^2}{\left(p+i\epsilon q\right)^2r^2c^2}}_{c,x}`$
$`={\displaystyle \frac{1}{4ϵ}}{\displaystyle \frac{\left(1\right)_{2ϵ}}{\left(\frac{1}{2}\right)_{2ϵ}}}F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,12ϵ,1}{22ϵ,1+ϵ}}|1\right)`$
$`{\displaystyle \frac{1}{2ϵ}}{\displaystyle \frac{\left(2\right)_{2ϵ}\left(1\right)_ϵ\left(\frac{1}{2}\right)_ϵ}{\left(1\right)_ϵ\left(\frac{1}{2}\right)_ϵ\left(\frac{3}{2}\right)_{2ϵ}}}F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{1}{2}ϵ,\frac{1}{2}2ϵ,1}{\frac{3}{2}2ϵ,\frac{1}{2}+ϵ}}|1\right)`$
$`+{\displaystyle \frac{1}{8ϵ\left(13ϵ\right)}}e^{i\pi ϵ}{\displaystyle \frac{\left(2\right)_{2ϵ}\left(1\right)_{2ϵ}\left(1\right)_ϵ\left(\frac{3}{2}\right)_ϵ}{\left(1\right)_ϵ\left(2\right)_{3ϵ}}}`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,13ϵ,\frac{3}{2}ϵ}{23ϵ,23ϵ}}|1\right).`$ (C.133)
Expanding in powers of $`ϵ`$ and then taking the real parts, we obtain

$`\mathrm{Re}s^2{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}{\displaystyle \frac{c^2}{\left(p+i\epsilon q\right)^2r^2c^2}}_{c,x}`$
$`={\displaystyle \frac{\pi ^2}{24}}+O\left(ϵ\right).`$ (C.134)
Inserting the sum of the integrals (C.131) and (C.134) into the thermal integral (C.113), we obtain

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{s^2n_F^2\left(s\right)}{p^2q^2}}\mathrm{Re}c^2{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}_c`$ $`=`$ $`O\left(ϵ\right).`$ (C.135)
It remains only to evaluate the integral in Eq. (C.124) with $`𝐩𝐪`$ in the numerator. We begin by using the identity

$`c^2{\displaystyle \frac{𝐩𝐪}{r^2}}{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}_{c,x}`$
$`={\displaystyle \frac{p^2+q^2}{\left(p^2q^2+i\epsilon \right)^2}}c^2_c{\displaystyle \frac{𝐩𝐪}{r^2}}_x`$
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{\pm }{}}{\displaystyle \frac{1}{\left(p+i\epsilon \pm q\right)^2}}{\displaystyle \frac{𝐩𝐪c^4}{\left(p+i\epsilon \pm q\right)^2r^2c^2}}_{c,x}.`$ (C.136)
In the first term on the right side, the average over $`c`$ is a simple multiplicative factor: $`c^2_c=1/\left(32ϵ\right)`$. The average over $`x`$ gives hypergeometric functions of argument $`\beta `$:

$$\frac{𝐩𝐪}{r^2}_x=\frac{1}{8}\beta \left[F\left(\genfrac{}{}{0pt}{}{1ϵ,1}{32ϵ}|\beta \right)F\left(\genfrac{}{}{0pt}{}{2ϵ,1}{32ϵ}|\beta \right)\right].$$
(C.137)
The integral over $`\beta `$ gives hypergeometric functions of argument 1:

$`s^2{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}{\displaystyle \frac{p^2+q^2}{\left(p^2q^2\right)^2}}{\displaystyle \frac{𝐩𝐪}{r^2}}_x`$
$`={\displaystyle \frac{1}{8}}{\displaystyle \frac{\left(2\right)_{2ϵ}}{\left(\frac{3}{2}\right)_{2ϵ}}}[F\left({\displaystyle \genfrac{}{}{0pt}{}{22ϵ,1ϵ,1}{\frac{3}{2}2ϵ,32ϵ}}\right|1)`$
$`F\left({\displaystyle \genfrac{}{}{0pt}{}{22ϵ,2ϵ,1}{\frac{3}{2}2ϵ,32ϵ}}\right|1)]+{\displaystyle \frac{1}{12}}{\displaystyle \frac{\left(3\right)_{2ϵ}}{\left(\frac{5}{2}\right)_{2ϵ}}}[F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,1}{\frac{5}{2}2ϵ}}\right|1)`$
$`F\left({\displaystyle \genfrac{}{}{0pt}{}{2ϵ,1}{\frac{5}{2}2ϵ}}\right|1)].`$ (C.138)
Expanding in powers of $`ϵ`$, we obtain

$$s^2_0^1𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}\frac{p^2+q^2}{\left(p^2q^2\right)^2}\frac{𝐩𝐪}{r^2}_x=\frac{\pi ^2}{16}+O\left(ϵ\right).$$
(C.139)
In the second term of (C.136), the average over $`c`$ is given by (C.127). In the $`+q`$ term, the average over $`x=\widehat{𝐩}\widehat{𝐪}`$ is

$`xF\left({\displaystyle \genfrac{}{}{0pt}{}{1,\frac{5}{2}}{\frac{7}{2}ϵ}}\right|{\displaystyle \frac{r^2}{\left(p+q\right)^2}})_x={\displaystyle \frac{52ϵ}{4ϵ}}[F\left({\displaystyle \genfrac{}{}{0pt}{}{2ϵ,1,\frac{5}{2}}{32ϵ,1+ϵ}}\right|\beta )`$
$`F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,1,\frac{5}{2}}{32ϵ,1+ϵ}}\right|\beta )]+{\displaystyle \frac{5}{4ϵ}}{\displaystyle \frac{\left(1\right)_ϵ\left(1\right)_{2ϵ}\left(3\right)_{2ϵ}\left(\frac{7}{2}\right)_ϵ}{\left(1\right)_ϵ\left(3\right)_{3ϵ}}}\beta ^ϵ`$
$`\times \left[F\left({\displaystyle \genfrac{}{}{0pt}{}{12ϵ,\frac{5}{2}ϵ}{33ϵ}}|\beta \right){\displaystyle \frac{12ϵ}{1ϵ}}F\left({\displaystyle \genfrac{}{}{0pt}{}{22ϵ,\frac{5}{2}ϵ}{33ϵ}}|\beta \right)\right].`$
Integrating over $`\beta `$, we obtain hypergeometric functions of argument 1:

$`{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}{\displaystyle \frac{𝐩𝐪c^4}{\left(p+q\right)^2r^2c^2}}_{c,x}`$
$`={\displaystyle \frac{1}{4ϵ\left(32ϵ\right)}}{\displaystyle \frac{\left(2\right)_{2ϵ}}{\left(\frac{5}{2}\right)_{2ϵ}}}[F\left({\displaystyle \genfrac{}{}{0pt}{}{22ϵ,2ϵ,1,\frac{5}{2}}{\frac{5}{2}2ϵ,32ϵ,1+ϵ}}\right|1)`$
$`F\left({\displaystyle \genfrac{}{}{0pt}{}{22ϵ,1ϵ,1,\frac{5}{2}}{\frac{5}{2}2ϵ,32ϵ,1+ϵ}}\right|1)]`$
$`+{\displaystyle \frac{1}{6ϵ\left(23ϵ\right)}}{\displaystyle \frac{\left(1\right)_ϵ\left(1\right)_{2ϵ}\left(3\right)_{2ϵ}\left(\frac{3}{2}\right)_ϵ}{\left(1\right)_ϵ\left(\frac{5}{2}\right)_{3ϵ}}}`$
$`\times [F\left({\displaystyle \genfrac{}{}{0pt}{}{23ϵ,12ϵ,\frac{5}{2}ϵ}{\frac{5}{2}3ϵ,33ϵ}}\right|1)`$
$`{\displaystyle \frac{12ϵ}{1ϵ}}F\left({\displaystyle \genfrac{}{}{0pt}{}{23ϵ,22ϵ,\frac{5}{2}ϵ}{\frac{5}{2}3ϵ,33ϵ}}\right|1)].`$ (C.141)
Expanding in powers of $`ϵ`$, we obtain

$`{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}{\displaystyle \frac{𝐩𝐪c^4}{\left(p+q\right)^2r^2c^2}}_{c,x}={\displaystyle \frac{\pi ^26}{18}}.`$
(C.142)
In the $`q`$ term in the integral of the second term of (C.136), we integrate over $`\beta `$ before averaging over $`x`$. The integral over $`\beta `$ can be expressed in terms of hypergeometric functions of type $`{}_{2}{}^{}F_{1}^{}`$:

$`s^2{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}{\displaystyle \frac{4𝐩𝐪}{\left(pq\right)^2}}{\displaystyle \frac{c^4}{\left(p+i\epsilon q\right)^2r^2c^2}}_c`$
$`={\displaystyle \frac{1}{2\left(32ϵ\right)ϵ}}{\displaystyle \frac{\left(2\right)_{2ϵ}}{\left(\frac{1}{2}\right)_{2ϵ}}}\left(12y\right)F\left({\displaystyle \genfrac{}{}{0pt}{}{22ϵ,1}{1+ϵ}}|1y\right)`$
$`{\displaystyle \frac{1}{4\left(32ϵ\right)ϵ}}{\displaystyle \frac{\left(1\right)_ϵ}{\left(\frac{1}{2}\right)_ϵ}}\left(12y\right)\left(1y\right)^{3/2}`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{1}{2}2ϵ,1}{\frac{1}{2}+ϵ}}|1y\right)`$
$`+{\displaystyle \frac{1}{8\left(23ϵ\right)ϵ}}e^{i\pi ϵ}\left(1\right)_ϵ\left(\frac{3}{2}\right)_ϵ\left(12y\right)\left(1y\right)^ϵ`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{23ϵ,\frac{5}{2}ϵ}{33ϵ}}|y\right).`$ (C.143)
The phase in the last term is $`e^{i\pi ϵ}`$ for the $`f(s_+,s_{},r)`$ term of (C.113), which comes from the $`p>q`$ region of the integral, and $`e^{i\pi ϵ}`$ for the $`f(s_{},s_+,r)`$ term, which comes from the $`p<q`$ region. The average over $`x=\widehat{𝐩}\widehat{𝐪}`$ can be expressed in terms of hypergeometric functions of type $`{}_{3}{}^{}F_{2}^{}`$ evaluated at 1:

$`s^2{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}{\displaystyle \frac{4𝐩𝐪}{\left(pq\right)^2}}{\displaystyle \frac{c^4}{\left(p+i\epsilon q\right)^2r^2c^2}}_{c,x}`$
$`={\displaystyle \frac{1}{4\left(32ϵ\right)ϵ}}{\displaystyle \frac{\left(2\right)_{2ϵ}}{\left(\frac{1}{2}\right)_{2ϵ}}}[F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,22ϵ,1}{32ϵ,1+ϵ}}\right|1)`$
$`F\left({\displaystyle \genfrac{}{}{0pt}{}{2ϵ,22ϵ,1}{32ϵ,1+ϵ}}\right|1)]`$
$`{\displaystyle \frac{1}{\left(32ϵ\right)ϵ}}{\displaystyle \frac{\left(1\right)_ϵ\left(3\right)_{2ϵ}\left(\frac{1}{2}\right)_ϵ}{\left(1\right)_ϵ\left(\frac{1}{2}\right)_ϵ\left(\frac{3}{2}\right)_{2ϵ}}}[F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{1}{2}ϵ,\frac{1}{2}2ϵ,1}{\frac{3}{2}2ϵ,\frac{1}{2}+ϵ}}\right|1)`$
$`+{\displaystyle \frac{1+2ϵ}{2\left(1ϵ\right)}}F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{1}{2}ϵ,\frac{1}{2}2ϵ,1}{\frac{3}{2}2ϵ,\frac{1}{2}+ϵ}}\right|1)]`$
$`+{\displaystyle \frac{1}{16\left(23ϵ\right)ϵ}}e^{i\pi ϵ}{\displaystyle \frac{\left(1\right)_ϵ\left(2\right)_{2ϵ}\left(2\right)_{2ϵ}\left(\frac{3}{2}\right)_ϵ}{\left(1\right)_ϵ\left(3\right)_{3ϵ}}}`$
$`\times [F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,23ϵ,\frac{5}{2}ϵ}{33ϵ,33ϵ}}\right|1)`$
$`{\displaystyle \frac{1ϵ}{12ϵ}}F\left({\displaystyle \genfrac{}{}{0pt}{}{2ϵ,23ϵ,\frac{5}{2}ϵ}{33ϵ,33ϵ}}\right|1)].`$ (C.144)
The expansion of the real part of the integral in powers of $`ϵ`$ is

$`s^2{\displaystyle _0^1}𝑑\beta \beta ^{2ϵ}\left(1\beta \right)^{1/2}`$
$`\times \mathrm{Re}{\displaystyle \frac{4𝐩𝐪}{\left(pq\right)^2}}{\displaystyle \frac{c^4}{\left(p+i\epsilon q\right)^2r^2c^2}}_{c,x}={\displaystyle \frac{9\pi ^2}{18}}+O\left(ϵ\right).`$
(C.145)
Inserting (C.139), (C.142), and (C.145) into the thermal integral of (C.136), we obtain

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{s^2n_F^2\left(s\right)}{p^2q^2}}{\displaystyle \frac{𝐩𝐪}{r^2}}\mathrm{Re}c^2{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}{\displaystyle \frac{\pi ^21}{6\pi ^2}}\left[{\displaystyle \frac{\pi ^2}{12}}\mathrm{log}2\right].`$ (C.146)
Inserting this along with (C.135) into (C.124), we obtain

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{s^2n_F^2\left(s\right)}{p^2r^2}}\mathrm{Re}c^2{\displaystyle \frac{r^2c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,rc)}}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}{\displaystyle \frac{1\pi ^2}{6\pi ^2}}\left[{\displaystyle \frac{\pi ^2}{12}}\mathrm{log}2\right].`$ (C.147)
Adding this integral to the subtracted integral in (C.111), we obtain the final result in (C.84). The subtracted integral appearing in (C.112) vanishes due to antisymmetry of the integrand. Thus the final result (C.85) is given by (C.112).

The integrals (C.86) and (C.87) can be computed directly in three dimensions, as described above. The integrals (C.88)– (C.91) are divergent and require subtractions to remove the divergences. We first isolate the divergent part which come from the region $`q0`$. We need one subtraction:

$`n_B\left(q\right)`$ $`=`$ $`\left(n_B\left(q\right){\displaystyle \frac{T}{q}}+{\displaystyle \frac{1}{2}}\right)+{\displaystyle \frac{T}{q}}{\displaystyle \frac{1}{2}}.`$ (C.148)
In the integral (C.89), it is convenient to first use the identity $`r_c^2=p^2+2𝐩𝐪/c+q^2/c^2`$ to expand it into 3 integrals, two of which are (C.86) and (C.88). In the third integral, the subtraction (C.148) is needed to remove the divergences.

For the convergent terms, the HTL average over $`c`$ and the angular average over $`x=\widehat{𝐩}\widehat{𝐪}`$ can be calculated in three dimensions:

$`\mathrm{Re}c{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}_{c,x}={\displaystyle \frac{1}{6\left(4p^2q^2\right)}}`$
$`+{\displaystyle \frac{q^2\left(4p^2+3q^2\right)}{3\left(4p^2q^2\right)^3}}\mathrm{log}{\displaystyle \frac{2p}{q}}+{\displaystyle \frac{\left(p+q\right)\left(4p^2+2pq+q^2\right)}{12pq\left(2p+q\right)^3}}\mathrm{log}{\displaystyle \frac{p+q}{p}}`$
$`{\displaystyle \frac{\left(pq\right)\left(4p^22pq+q^2\right)}{12pq\left(2pq\right)^3}}\mathrm{log}{\displaystyle \frac{\left|pq\right|}{p}},`$ (C.149)
$`\mathrm{Re}\widehat{𝐩}\widehat{𝐪}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}_{c,x}={\displaystyle \frac{1}{6pq}}{\displaystyle \frac{q\left(12p^2q^2\right)}{6p\left(4p^2q^2\right)^2}}\mathrm{log}{\displaystyle \frac{4p}{q}}`$
$`+{\displaystyle \frac{\left(p+q\right)\left(2p^22pqq^2\right)}{12p^2q\left(2p+q\right)^2}}\mathrm{log}{\displaystyle \frac{p+q}{4p}}`$
$`+{\displaystyle \frac{\left(pq\right)\left(2p^2+2pqq^2\right)}{12p^2q\left(2pq\right)^2}}\mathrm{log}{\displaystyle \frac{\left|pq\right|}{4p}},`$ (C.150)
$`\text{Re}{\displaystyle \frac{r_c^2p^2}{q^2}}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+iϵ,q,r_c)}}c^1{\displaystyle \frac{1}{q^2}}c^1+{\displaystyle \frac{\mathrm{log}2}{q^2}}_{c,x}`$
$`={\displaystyle \frac{1}{4pq^2}}\left[q\mathrm{log}{\displaystyle \frac{p+q}{\left|pq\right|}}+p\mathrm{log}{\displaystyle \frac{\left|p^2q^2\right|}{p^2}}\right].`$ (C.151)
The remaining 2-dimensional integral over $`p`$ and $`q`$ can then be evaluated numerically:

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}\left({\displaystyle \frac{n_B\left(q\right)}{q}}{\displaystyle \frac{T}{q^2}}+{\displaystyle \frac{1}{2q}}\right){\displaystyle \frac{p^2}{q^2}}`$
$`\times \mathrm{Re}c{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}[1.480\times 10^2],`$
(C.152)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}\left({\displaystyle \frac{n_B\left(q\right)}{q}}{\displaystyle \frac{T}{q^2}}+{\displaystyle \frac{1}{2q}}\right){\displaystyle \frac{𝐩𝐪}{q^2}}`$
$`\times \mathrm{Re}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}[2.832\times 10^3],`$
(C.153)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}\text{Re}{\displaystyle \frac{r_c^2p^2}{q^2}}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+iϵ,q,r_c)}}c^1`$
$`{\displaystyle \frac{1}{q^2}}c^1+{\displaystyle \frac{\mathrm{log}2}{q^2}}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}[4.134\times 10^2],`$ (C.154)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_B\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}\text{Re}{\displaystyle \frac{r_c^2p^2}{q^2}}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+iϵ,q,r_c)}}c^1`$
$`{\displaystyle \frac{1}{q^2}}c^1+{\displaystyle \frac{\mathrm{log}2}{q^2}}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}[2.530\times 10^1].`$ (C.155)
The integrals involving the terms subtracted from $`n\left(q\right)`$ in (C.148) are divergent, so the HTL average over $`c`$ and the angular average over $`x=\widehat{𝐩}\widehat{𝐪}`$ must be calculated in $`32ϵ`$ dimensions. The first step in the calculation of the subtracted terms is to replace the average over $`c`$ of the integral over $`q`$ by an average over $`c`$ and $`x`$:

$`{\displaystyle _𝐪}{\displaystyle \frac{1}{q^n}}f\left(c\right){\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}_c`$
$`=\left(1\right)^{n1}{\displaystyle \frac{1}{8\pi ^2ϵ}}{\displaystyle \frac{\left(1\right)_{2ϵ}\left(1\right)_{2ϵ}}{\left(\frac{3}{2}\right)_ϵ}}\left(e^\gamma \mu ^2\right)^ϵ\left(2p\right)^{1n2ϵ}`$
$`\times f\left(c\right)c^{3n2ϵ}\left(1c^2\right)^{n2+2ϵ}{\displaystyle \underset{\pm }{}}\left(xci\epsilon \right)^{1n2ϵ}_{c,x}.`$
(C.156)
The integral over $`p`$ can now be evaluated easily using either (B.54) or

$$_𝐩n_F\left(p\right)p^{22ϵ}=\frac{1}{2\pi ^2}\frac{\left(1\right)_{4ϵ}}{\left(\frac{3}{2}\right)_ϵ}\left(12^{4ϵ}\right)\zeta \left(14ϵ\right)\left(e^\gamma \mu ^2\right)^ϵT^{14ϵ}.$$
(C.157)
It remains only to calculate the averages over $`c`$ and $`x`$. The averages over $`x`$ give $`{}_{2}{}^{}F_{1}^{}`$ hypergeometric functions with argument $`\left[\left(1c\right)/2i\epsilon \right]^1`$:

$`\left(xci\epsilon \right)^{n2ϵ}_x=\left(1c\right)^{n2ϵ}`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,n+2ϵ}{22ϵ}}|\left[\left(1c\right)/2i\epsilon \right]^1\right),`$ (C.158)
$`x\left(xci\epsilon \right)^{n2ϵ}_x={\displaystyle \frac{1}{2}}\left(1c\right)^{n2ϵ}`$
$`\times [F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,n+2ϵ}{32ϵ}}\right|[(1c)/2i\epsilon ]^1)`$
$`F\left({\displaystyle \genfrac{}{}{0pt}{}{2ϵ,n+2ϵ}{32ϵ}}\right|[(1c)/2i\epsilon ]^1)].`$ (C.159)
Using a transformation formula, the arguments can be changed to $`\left(1c\right)/2i\epsilon `$. If the expressions (C.158) and (C.159) are averaged over $`c`$ with a weight that is an even function of $`c`$, the $`+`$ and $``$ terms combine to give $`{}_{3}{}^{}F_{2}^{}`$ hypergeometric functions with argument 1. For example,

$`\left(1c^2\right)^{2ϵ}{\displaystyle \underset{\pm }{}}\left(xci\epsilon \right)^{12ϵ}_{c,x}`$
$`={\displaystyle \frac{1}{3ϵ}}{\displaystyle \frac{\left(2\right)_{2ϵ}\left(1\right)_ϵ\left(\frac{3}{2}\right)_ϵ}{\left(1\right)_ϵ\left(1\right)_ϵ}}`$
$`\times \{e^{i\pi ϵ}{\displaystyle \frac{\left(1\right)_{3ϵ}\left(1\right)_{2ϵ}}{\left(1\right)_{2ϵ}\left(2\right)_ϵ}}F\left({\displaystyle \genfrac{}{}{0pt}{}{12ϵ,1ϵ,ϵ}{2ϵ,13ϵ}}\right|1)`$
$`+e^{i2\pi ϵ}{\displaystyle \frac{\left(1\right)_{3ϵ}\left(1\right)_ϵ}{\left(1\right)_{4ϵ}\left(2\right)_{2ϵ}}}F\left({\displaystyle \genfrac{}{}{0pt}{}{1+ϵ,1+2ϵ,4ϵ}{2+2ϵ,1+3ϵ}}\right|1)\}.`$ (C.160)
Upon expanding the hypergeometric functions in powers of $`ϵ`$ and taking the real parts, we obtain

$`\mathrm{Re}\left(1c^2\right)^{2ϵ}{\displaystyle \underset{\pm }{}}\left(xci\epsilon \right)^{12ϵ}_{c,x}`$
$`=\pi ^2\left[ϵ+2\left(1\mathrm{log}2\right)ϵ^2\right],`$ (C.161)
$`\mathrm{Re}c^2\left(1c^2\right)^{2ϵ}{\displaystyle \underset{\pm }{}}\left(xci\epsilon \right)^{12ϵ}_{c,x}`$
$`=\pi ^2\left[{\displaystyle \frac{1}{3}}ϵ+{\displaystyle \frac{2}{9}}\left(23\mathrm{log}2\right)ϵ^2\right],`$ (C.162)
$`\mathrm{Re}\left(1c^2\right)^{2+2ϵ}{\displaystyle \underset{\pm }{}}\left(xci\epsilon \right)^{32ϵ}_{c,x}`$
$`=\pi ^2\left[{\displaystyle \frac{8}{3}}ϵ^2\right],`$ (C.163)
$`\mathrm{Re}x\left(1c^2\right)^{1+2ϵ}{\displaystyle \underset{\pm }{}}\left(xci\epsilon \right)^{22ϵ}_{c,x}`$
$`=\pi ^2\left[{\displaystyle \frac{2}{3}}ϵ+{\displaystyle \frac{2}{9}}\left(16\mathrm{log}2\right)ϵ^2\right].`$ (C.164)
If the expressions (C.158) and (C.159) are averaged over $`c`$ with a weight that is an odd function of $`c`$, they reduce to integrals of $`{}_{2}{}^{}F_{1}^{}`$ hypergeometric functions with argument $`y`$. For example,

$`c\left(1c^2\right)^{1+2ϵ}{\displaystyle \underset{\pm }{}}\left(xci\epsilon \right)^{22ϵ}_{c,x}={\displaystyle \frac{\left(2\right)_{2ϵ}\left(\frac{3}{2}\right)_ϵ}{\left(1\right)_ϵ\left(1\right)_ϵ}}`$
$`\times \{2e^{i\pi ϵ}{\displaystyle \frac{\left(1\right)_{3ϵ}}{\left(2\right)_{2ϵ}}}{\displaystyle _0^1}dyy^{2ϵ}(1y)^{1+ϵ}|12y|`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,ϵ}{3ϵ}}|y\right)`$
$`{\displaystyle \frac{8}{3\left(1+3ϵ\right)}}e^{2i\pi ϵ}{\displaystyle \frac{\left(1\right)_{3ϵ}}{\left(1\right)_{4ϵ}}}{\displaystyle _0^1}𝑑yy^{1+ϵ}\left(1y\right)^{1+ϵ}\left|12y\right|`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{2+2ϵ,1+4ϵ}{2+3ϵ}}\right|y)\}.`$ (C.165)
The resulting expansions for the real parts of the averages over $`c`$ and $`x`$ are

$`\mathrm{Re}c\left(1c^2\right)^{1+2ϵ}{\displaystyle \underset{\pm }{}}\left(xci\epsilon \right)^{22ϵ}_{c,x}`$
$`=1+{\displaystyle \frac{14\left(1\mathrm{log}2\right)}{3}}ϵ,`$ (C.166)
$`\mathrm{Re}xc\left(1c^2\right)^{2ϵ}{\displaystyle \underset{\pm }{}}\left(xci\epsilon \right)^{12ϵ}_{c,x}={\displaystyle \frac{2\left(1\mathrm{log}2\right)}{3}}`$
$`+\left({\displaystyle \frac{4}{9}}+{\displaystyle \frac{8}{9}}\mathrm{log}2{\displaystyle \frac{4}{3}}\mathrm{log}^22+{\displaystyle \frac{\pi ^2}{18}}\right)ϵ.`$ (C.167)
Multiplying each of these expansions by the appropriate factors from the integral over $`q`$ in (C.156) and the integral over $`p`$ in (C.157) or (B.54), we obtain

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{p^2}{q^3}}\mathrm{Re}c^{1+2ϵ}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$
$`\times \left({\displaystyle \frac{1}{48}}\right)\left[{\displaystyle \frac{1}{ϵ}}{\displaystyle \frac{2}{3}}{\displaystyle \frac{4}{3}}\mathrm{log}2+4{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right],`$ (C.168)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{p^2}{q^4}}\mathrm{Re}c^{1+2ϵ}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}_c=O\left(ϵ\right),`$ (C.169)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{𝐩𝐪}{q^3}}\mathrm{Re}c^{2ϵ}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$
$`\times ({\displaystyle \frac{1}{36}})[(1\mathrm{log}2)({\displaystyle \frac{1}{ϵ}}+{\displaystyle \frac{14}{3}}4\mathrm{log}2`$
$`+4{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}})+{\displaystyle \frac{\pi ^2}{12}}],`$ (C.170)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{𝐩𝐪}{q^4}}\mathrm{Re}c^{2ϵ}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}_c`$
$`={\displaystyle \frac{T}{\left(4\pi \right)^2}}\left({\displaystyle \frac{1}{6}}\mathrm{log}2\right).`$ (C.171)
Adding Eq. (C.169) to the subtracted integral (C.152) we obtain the final result in Eq. (C.88). Combining (C.153) with (C.170) and (C.171), we obtain

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)n_B\left(q\right)}{pq}}{\displaystyle \frac{𝐩𝐪}{q^2}}\mathrm{Re}c^{2ϵ}{\displaystyle \frac{r_c^2p^2q^2}{\mathrm{\Delta }(p+i\epsilon ,q,r_c)}}_c`$
$`={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}\left({\displaystyle \frac{1\mathrm{log}2}{72}}\right)\left[{\displaystyle \frac{1}{ϵ}}15.2566\right].`$ (C.172)
The integral (C.89) is obtained from (C.86), (C.88) and (C.172). Finally consider (C.90) and (C.91). In order to evaluate them we need two subtractions for each integral

$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(p\right)}{q}}{\displaystyle \frac{1}{q^2}}c^{2ϵ}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$
$`\times \left({\displaystyle \frac{1}{12}}\right)\left[{\displaystyle \frac{1}{ϵ}}+2+2\mathrm{log}2+2\gamma +2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right],`$ (C.173)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_F\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{1}{q^2}}c^{1+2ϵ}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$
$`\times ({\displaystyle \frac{1}{24}})[{\displaystyle \frac{1}{ϵ^2}}+(2+2\gamma +4\mathrm{log}2+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}){\displaystyle \frac{1}{ϵ}}`$
$`+53.1064],`$ (C.174)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_B\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{1}{q^2}}c^{2ϵ}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$
$`\times \left({\displaystyle \frac{1}{6}}\right)\left[{\displaystyle \frac{1}{ϵ}}+2+4\mathrm{log}2+2\gamma +2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}\right],`$ (C.175)
$`{\displaystyle _{\mathrm{𝐩𝐪}}}{\displaystyle \frac{n_B\left(p\right)}{p}}{\displaystyle \frac{n_F\left(q\right)}{q}}{\displaystyle \frac{1}{q^2}}c^{1+2ϵ}_c={\displaystyle \frac{T^2}{\left(4\pi \right)^2}}\left({\displaystyle \frac{\mu }{4\pi T}}\right)^{4ϵ}`$
$`\times ({\displaystyle \frac{1}{12}})[{\displaystyle \frac{1}{ϵ^2}}+(2+2\gamma +6\mathrm{log}2+2{\displaystyle \frac{\zeta ^{}\left(1\right)}{\zeta \left(1\right)}}){\displaystyle \frac{1}{ϵ}}`$
$`+69.7097].`$ (C.176)
The subtractions can be evaluated directly in three dimensions and the results are given in Eqs. (C.154)–(C.155) The integrals (C.90) and (C.91) are then given by the by the sum of the difference terms (C.154) and (C.155) and the subtraction terms (C.173)–(C.176).

### C.3 4-dimensional integrals

In the sum-integral formula (B.39), the second term on the right side involves an integral over 4-dimensional Euclidean momenta. The integrands are functions of the integration variable $`Q`$ and $`R=\left(P+Q\right)`$. The simplest integrals to evaluate are those whose integrands are independent of $`P_0`$:

$`{\displaystyle _Q}{\displaystyle \frac{1}{Q^2r^2}}`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{2ϵ}\mathrm{\hspace{0.33em}2}\left[{\displaystyle \frac{1}{ϵ}}+42\mathrm{log}2\right],`$ (C.177)
$`{\displaystyle _Q}{\displaystyle \frac{q^2}{Q^2r^4}}`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{2ϵ}\mathrm{\hspace{0.33em}2}\left[{\displaystyle \frac{1}{ϵ}}+12\mathrm{log}2\right],`$ (C.178)
$`{\displaystyle _Q}{\displaystyle \frac{1}{Q^2r^4}}`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{22ϵ}\left(2\right)\left[1+\left(22\mathrm{log}2\right)ϵ\right].`$
Another simple integral that is needed depends only on $`P^2=P_0^2+p^2`$:

$`{\displaystyle _Q}{\displaystyle \frac{1}{Q^2R^2}}={\displaystyle \frac{1}{\left(4\pi \right)^2}}\left(e^\gamma \mu ^2\right)^ϵ\left(P^2\right)^ϵ{\displaystyle \frac{1}{ϵ}}{\displaystyle \frac{\left(1\right)_ϵ\left(1\right)_ϵ\left(1\right)_ϵ}{\left(2\right)_{2ϵ}}},`$ (C.180)
where $`\left(a\right)_b`$ is Pochhammer’s symbol which is defined in (C.212). We need the following weighted averages over $`c`$ of this function evaluated at $`P=(ip,𝐩/c)`$:

$`c^{1+2ϵ}{\displaystyle _Q}{\displaystyle \frac{1}{Q^2R^2}}|_{P(ip,𝐩/c)}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{2ϵ}`$
$`\times {\displaystyle \frac{1}{4}}\left[{\displaystyle \frac{1}{ϵ^2}}+{\displaystyle \frac{2\mathrm{log}2}{ϵ}}+2\mathrm{log}^22+{\displaystyle \frac{3\pi ^2}{4}}\right],`$ (C.181)
$`c^{1+2ϵ}{\displaystyle _Q}{\displaystyle \frac{1}{Q^2R^2}}|_{P(ip,𝐩/c)}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{2ϵ}`$
$`\times {\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{1}{ϵ}}+2\mathrm{log}2\right].`$ (C.182)
The remaining integrals are functions of $`P_0`$ that must be analytically continued to the point $`P_0=ip+\epsilon `$. Several of these integrals are straightforward to evaluate:

$`{\displaystyle _Q}{\displaystyle \frac{q^2}{Q^2R^2}}|_{P_0=ip}=0,`$ (C.183)
$`{\displaystyle _Q}{\displaystyle \frac{q^2}{Q^2r^2R^2}}|_{P_0=ip}={\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{2ϵ}`$
$`\times (1)[{\displaystyle \frac{1}{ϵ^2}}+{\displaystyle \frac{12\mathrm{log}2}{ϵ}}+102\mathrm{log}2`$
$`+2\mathrm{log}^22{\displaystyle \frac{7\pi ^2}{12}}],`$ (C.184)
$`{\displaystyle _Q}{\displaystyle \frac{1}{Q^2r^2R^2}}|_{P_0=ip}={\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{22ϵ}\left[{\displaystyle \frac{1}{ϵ}}22\mathrm{log}2\right].`$
We also need a weighted average over $`c`$ of the integral in (C.183) evaluated at $`P=(ip,𝐩/c)`$. The integral itself is

$`{\displaystyle _Q}{\displaystyle \frac{q^2}{Q^2R^2}}|_{P(ip,𝐩/c)}={\displaystyle \frac{1}{\left(4\pi \right)^2}}\left(e^\gamma \mu ^2\right)^ϵp^{22ϵ}{\displaystyle \frac{\left(1\right)_ϵ}{ϵ}}`$
$`\times {\displaystyle \frac{1}{4}}{\displaystyle \frac{\left(1\right)_ϵ\left(1\right)_ϵ}{\left(2\right)_{2ϵ}}}\left({\displaystyle \frac{1}{32ϵ}}+c^2\right)c^{2+2ϵ}\left(1c^2\right)^ϵ.`$ (C.186)
The weighted averages are

$`c^{1+2ϵ}{\displaystyle _Q}{\displaystyle \frac{q^2}{Q^2R^2}}|_{P(ip,𝐩/c)}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{22ϵ}`$
$`\times {\displaystyle \frac{1}{48}}[{\displaystyle \frac{1}{ϵ^2}}+{\displaystyle \frac{2\left(10+3\mathrm{log}2\right)}{3ϵ}}`$
$`+{\displaystyle \frac{4}{9}}+{\displaystyle \frac{40}{3}}\mathrm{log}2+2\mathrm{log}^22+{\displaystyle \frac{3\pi ^2}{4}}].`$ (C.187)
$`c^{1+2ϵ}{\displaystyle _Q}{\displaystyle \frac{q^2}{Q^2R^2}}|_{P(ip,𝐩/c)}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{22ϵ}`$
$`\times {\displaystyle \frac{1}{16}}\left[{\displaystyle \frac{1}{ϵ^2}}+{\displaystyle \frac{2\mathrm{log}2}{ϵ}}+2\mathrm{log}^22+{\displaystyle \frac{3\pi ^2}{4}}\right].`$ (C.188)
The most difficult 4-dimensional integrals to evaluate involve an HTL average of an integral with denominator $`R_0^2+r^2c^2`$:

$`\mathrm{Re}{\displaystyle _Q}{\displaystyle \frac{1}{Q^2}}{\displaystyle \frac{c^2}{R_0^2+r^2c^2}}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{2ϵ}`$
$`\times \left[{\displaystyle \frac{22\mathrm{log}2}{ϵ}}+84\mathrm{log}2+4\mathrm{log}^22{\displaystyle \frac{\pi ^2}{2}}\right],`$ (C.189)
$`\mathrm{Re}{\displaystyle _Q}{\displaystyle \frac{1}{Q^2}}{\displaystyle \frac{c^2\left(1c^2\right)}{R_0^2+r^2c^2}}_c{\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{2ϵ}`$
$`={\displaystyle \frac{1}{3}}\left[{\displaystyle \frac{1}{ϵ}}+{\displaystyle \frac{20}{3}}6\mathrm{log}2\right],`$ (C.190)
$`\mathrm{Re}{\displaystyle _Q}{\displaystyle \frac{1}{Q^2}}{\displaystyle \frac{c^4}{R_0^2+r^2c^2}}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{2ϵ}`$
$`\times \left[{\displaystyle \frac{56\mathrm{log}2}{3ϵ}}+{\displaystyle \frac{52}{9}}2\mathrm{log}2+4\mathrm{log}^22{\displaystyle \frac{\pi ^2}{2}}\right],`$ (C.191)
$`\mathrm{Re}{\displaystyle _Q}{\displaystyle \frac{1}{Q^2r^2}}{\displaystyle \frac{c^2}{R_0^2+r^2c^2}}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{22ϵ}`$
$`\times \left({\displaystyle \frac{1}{4}}\right)\left[{\displaystyle \frac{1}{ϵ}}+{\displaystyle \frac{4}{3}}+{\displaystyle \frac{2}{3}}\mathrm{log}2\right],`$ (C.192)
$`\mathrm{Re}{\displaystyle _Q}{\displaystyle \frac{q^2}{Q^2r^2}}{\displaystyle \frac{c^2}{R_0^2+r^2c^2}}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{2ϵ}`$
$`\times \left[{\displaystyle \frac{1316\mathrm{log}2}{12ϵ}}+{\displaystyle \frac{29}{9}}{\displaystyle \frac{19}{18}}\mathrm{log}2+{\displaystyle \frac{8}{3}}\mathrm{log}^22{\displaystyle \frac{4}{9}}\pi ^2\right],`$
(C.193)
$`{\displaystyle _Q}{\displaystyle \frac{q^2p^2}{Q^2r^2\left(R_0^2+r^2c^2\right)}}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{2ϵ}\left[{\displaystyle \frac{\pi ^2}{3}}\right].`$ (C.194)
The analytic continuation to $`P_0=ip+\epsilon `$ is implied in these integrals and in all the 4-dimensional integrals in the remainder of this subsection.

We proceed to describe the evaluation of the integrals (C.189) and (C.191). The integral over $`Q_0`$ can be evaluated by introducing a Feynman parameter to combine $`Q^2`$ and $`R_0^2+r^2c^2`$ into a single denominator:

$`{\displaystyle _Q}{\displaystyle \frac{1}{Q^2\left(R_0^2+r^2c^2\right)}}={\displaystyle \frac{1}{4}}{\displaystyle _0^1}dx{\displaystyle _𝐫}[(1x+xc^2)r^2`$
$`+2(1x)𝐫𝐩+(1x)^2p^2i\epsilon ]^{3/2},`$ (C.195)
where we have carried out the analytic continuation to $`P_0=ip+\epsilon `$. Integrating over $`𝐫`$ and then over the Feynman parameter, we get a $`{}_{2}{}^{}F_{1}^{}`$ hypergeometric function with argument $`1c^2`$:

$`{\displaystyle _Q}{\displaystyle \frac{1}{Q^2\left(R_0^2+r^2c^2\right)}}={\displaystyle \frac{1}{\left(4\pi \right)^2}}\left(e^\gamma \mu ^2\right)^ϵp^{2ϵ}{\displaystyle \frac{\left(1\right)_ϵ}{ϵ}}`$
$`\times e^{i\pi ϵ}{\displaystyle \frac{\left(1\right)_{2ϵ}\left(1\right)_ϵ}{\left(2\right)_{3ϵ}}}\left(1c^2\right)^ϵF\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{3}{2}2ϵ,1ϵ}{23ϵ}}|1c^2\right).`$
(C.196)
The subsequent weighted averages over $`c`$ give $`{}_{3}{}^{}F_{2}^{}`$ hypergeometric functions with argument $`1`$:

$`{\displaystyle _Q}{\displaystyle \frac{1}{Q^2}}{\displaystyle \frac{c^2}{R_0^2+r^2c^2}}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\left(e^\gamma \mu ^2\right)^ϵp^{2ϵ}{\displaystyle \frac{\left(1\right)_ϵ}{ϵ}}`$
$`\times {\displaystyle \frac{1}{3}}e^{i\pi ϵ}{\displaystyle \frac{\left(\frac{3}{2}\right)_ϵ\left(1\right)_{2ϵ}\left(1\right)_{2ϵ}}{\left(\frac{5}{2}\right)_{2ϵ}\left(2\right)_{3ϵ}}}F\left({\displaystyle \genfrac{}{}{0pt}{}{12ϵ,\frac{3}{2}2ϵ,1ϵ}{\frac{5}{2}2ϵ,23ϵ}}|1\right),`$
(C.197)
$`{\displaystyle _Q}{\displaystyle \frac{1}{Q^2}}{\displaystyle \frac{c^2\left(1c^2\right)}{R_0^2+r^2c^2}}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\left(e^\gamma \mu ^2\right)^ϵp^{2ϵ}{\displaystyle \frac{\left(1\right)_ϵ}{ϵ}}`$
$`\times {\displaystyle \frac{2}{15}}e^{i\pi ϵ}{\displaystyle \frac{\left(\frac{3}{2}\right)_ϵ\left(1\right)_{2ϵ}\left(2\right)_{2ϵ}}{\left(\frac{7}{2}\right)_{2ϵ}\left(2\right)_{3ϵ}}}F\left({\displaystyle \genfrac{}{}{0pt}{}{22ϵ\frac{3}{2}2ϵ,1ϵ}{\frac{7}{2}2ϵ,23ϵ}}|1\right).`$
After expanding in powers of $`ϵ`$, the real part is (C.191).

The integral (C.192) has a factor of $`1/r^2`$ in the integrand. After using (C.195), it is convenient to use a second Feynman parameter to combine $`\left(1x+xc^2\right)r^2`$ with the other denominator before integrating over $`𝐫`$:

$`{\displaystyle _Q}{\displaystyle \frac{1}{Q^2r^2\left(R_0^2+r^2c^2\right)}}={\displaystyle \frac{3}{8}}{\displaystyle _0^1}𝑑x\left(1x+xc^2\right){\displaystyle _0^1}𝑑yy^{1/2}`$
$`\times {\displaystyle _𝐫}[(1x+xc^2)r^2+2y(1x)𝐫𝐩`$
$`+y(1x)^2p^2i\epsilon ]^{5/2}.`$ (C.199)
After integrating over $`𝐫`$ and then $`y`$, we obtain $`{}_{2}{}^{}F_{1}^{}`$ hypergeometric functions with arguments $`x\left(1c^2\right)`$. The integral over $`x`$ gives a $`{}_{2}{}^{}F_{1}^{}`$ hypergeometric function with argument $`1c^2`$:

$`{\displaystyle _Q}{\displaystyle \frac{1}{Q^2r^2\left(R_0^2+r^2c^2\right)}}={\displaystyle \frac{1}{\left(4\pi \right)^2}}\left(e^\gamma \mu ^2\right)^ϵp^{22ϵ}{\displaystyle \frac{\left(1\right)_ϵ}{ϵ}}`$
$`\times \{{\displaystyle \frac{\left(\frac{1}{2}\right)_ϵ\left(1\right)_ϵ}{\left(\frac{1}{2}\right)_{2ϵ}}}{\displaystyle \frac{3}{2\left(1+2ϵ\right)}}e^{i\pi ϵ}{\displaystyle \frac{\left(1\right)_{2ϵ}\left(1\right)_ϵ}{\left(1\right)_{3ϵ}}}(1c^2)^ϵ`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{1}{2}2ϵ,ϵ}{3ϵ}}\right|1c^2)\}.`$ (C.200)
After averaging over $`c`$, we get a hypergeometric functions with argument 1:

$`{\displaystyle _Q}{\displaystyle \frac{1}{Q^2r^2}}{\displaystyle \frac{c^2}{R_0^2+r^2c^2}}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\left(e^\gamma \mu ^2\right)^ϵp^{22ϵ}{\displaystyle \frac{\left(1\right)_ϵ}{ϵ}}`$
$`\times \{{\displaystyle \frac{1}{32ϵ}}{\displaystyle \frac{\left(\frac{1}{2}\right)_ϵ\left(1\right)_ϵ}{\left(\frac{1}{2}\right)_{2ϵ}}}{\displaystyle \frac{1}{2}}e^{i\pi ϵ}{\displaystyle \frac{\left(\frac{1}{2}\right)_ϵ\left(1\right)_{2ϵ}\left(2\right)_{2ϵ}}{\left(\frac{5}{2}\right)_{2ϵ}\left(1\right)_{3ϵ}}}`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{12ϵ,\frac{1}{2}2ϵ,ϵ}{\frac{5}{2}2ϵ,3ϵ}}\right|1)\}.`$ (C.201)
$`{\displaystyle _Q}{\displaystyle \frac{1}{Q^2r^2}}{\displaystyle \frac{1}{R_0^2+r^2c^2}}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\left(e^\gamma \mu ^2\right)^ϵp^{22ϵ}{\displaystyle \frac{\left(1\right)_ϵ}{ϵ}}`$
$`\times \{{\displaystyle \frac{\left(\frac{1}{2}\right)_ϵ\left(\frac{3}{2}\right)_ϵ}{\left(\frac{1}{2}\right)_{2ϵ}}}{\displaystyle \frac{1}{2}}e^{i\pi ϵ}{\displaystyle \frac{\left(1\right)_{2ϵ}^2\left(1\right)_ϵ}{\left(\frac{3}{2}\right)_{2ϵ}\left(1\right)_{3ϵ}}}`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{12ϵ,\frac{1}{2}2ϵ,ϵ}{\frac{3}{2}2ϵ,3ϵ}}\right|1)\}.`$ (C.202)
After expanding in powers of $`ϵ`$, the real part is (C.192).

To evaluate the integral (C.193), it is convenient to first express it as the sum of 3 integrals by expanding the factor of $`q^2`$ in the numerator as $`q^2=p^2+2𝐩𝐫+r^2`$:

$`{\displaystyle _Q}{\displaystyle \frac{q^2}{Q^2r^2\left(R_0^2+r^2c^2\right)}}={\displaystyle _Q}\left({\displaystyle \frac{p^2}{r^2}}+2{\displaystyle \frac{𝐩𝐫}{r^2}}+1\right)`$
$`\times {\displaystyle \frac{1}{Q^2\left(R_0^2+r^2c^2\right)}}.`$ (C.203)
To evaluate the integral with $`𝐩𝐫`$ in the numerator, we first combine the denominators using Feynman parameters as in (C.199). After integrating over $`𝐫`$ and then $`y`$, we obtain $`{}_{2}{}^{}F_{1}^{}`$ hypergeometric functions with arguments $`x\left(1c^2\right)`$. The integral over $`x`$ gives $`{}_{2}{}^{}F_{1}^{}`$ hypergeometric functions with arguments $`1c^2`$:

$`{\displaystyle _Q}{\displaystyle \frac{𝐩𝐫}{Q^2r^2\left(R_0^2+r^2c^2\right)}}={\displaystyle \frac{1}{\left(4\pi \right)^2}}\left(e^\gamma \mu ^2\right)^ϵp^{2ϵ}{\displaystyle \frac{\left(1\right)_ϵ}{2ϵ^2}}`$
$`\times \{{\displaystyle \frac{\left(\frac{3}{2}\right)_ϵ\left(1\right)_ϵ}{\left(\frac{3}{2}\right)_{2ϵ}}}+e^{i\pi ϵ}{\displaystyle \frac{\left(1\right)_{2ϵ}\left(1\right)_ϵ}{\left(1\right)_{3ϵ}}}(1c^2)^ϵ`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{\frac{3}{2}2ϵ,ϵ}{13ϵ}}\right|1c^2)\}.`$ (C.204)
After averaging over $`c`$, we get a hypergeometric function with argument 1:

$`{\displaystyle _Q}{\displaystyle \frac{𝐩𝐫}{Q^2r^2}}{\displaystyle \frac{c^2}{R_0^2+r^2c^2}}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\left(e^\gamma \mu ^2\right)^ϵp^{2ϵ}{\displaystyle \frac{\left(1\right)_ϵ}{2ϵ^2}}`$
$`\times \{{\displaystyle \frac{1}{32ϵ}}{\displaystyle \frac{\left(\frac{3}{2}\right)_ϵ\left(1\right)_ϵ}{\left(\frac{3}{2}\right)_{2ϵ}}}+{\displaystyle \frac{1}{3}}e^{i\pi ϵ}{\displaystyle \frac{\left(\frac{3}{2}\right)_ϵ\left(1\right)_{2ϵ}\left(1\right)_{2ϵ}}{\left(\frac{5}{2}\right)_{2ϵ}\left(1\right)_{3ϵ}}}`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{12ϵ,\frac{3}{2}2ϵ,ϵ}{\frac{5}{2}2ϵ,13ϵ}}\right|1)\}.`$ (C.205)
After expanding in powers of $`ϵ`$, the real part is

$`\mathrm{Re}{\displaystyle _Q}{\displaystyle \frac{𝐩𝐫}{Q^2r^2}}{\displaystyle \frac{c^2}{R_0^2+r^2c^2}}_c={\displaystyle \frac{1}{\left(4\pi \right)^2}}\mu ^{2ϵ}p^{2ϵ}[{\displaystyle \frac{1+\mathrm{log}2}{3ϵ}}`$
$`{\displaystyle \frac{20}{9}}+{\displaystyle \frac{14}{9}}\mathrm{log}2{\displaystyle \frac{2}{3}}\mathrm{log}^22+{\displaystyle \frac{\pi ^2}{36}}].`$ (C.206)
Combining this with (C.189) and (C.191), we obtain the integral (C.193).

To evaluate the integral (C.194), we first express the numerator as a sum of two integrals whose averages have been calculated:

$`{\displaystyle _Q}{\displaystyle \frac{q^2p^2}{Q^2r^2\left(R_0^2+r^2c^2\right)}}_x={\displaystyle _Q}{\displaystyle \frac{2𝐩𝐫+r^2}{Q^2r^2\left(R_0^2+r^2c^2\right)}}_x`$
$`={\displaystyle \frac{1}{\left(4\pi \right)^2}}\left(e^\gamma \mu ^2\right)^ϵp^{2ϵ}{\displaystyle \frac{\left(1\right)_ϵ}{ϵ}}\{{\displaystyle \frac{1}{ϵ}}{\displaystyle \frac{\left(\frac{3}{2}\right)_ϵ\left(1\right)_ϵ}{\left(\frac{3}{2}\right)_{2ϵ}}}`$
$`+e^{i\pi ϵ}{\displaystyle \frac{\left(1\right)_ϵ\left(1\right)_{2ϵ}}{\left(1\right)_{3ϵ}}}{\displaystyle \frac{1}{ϵ}}\left(1c^2\right)^ϵF\left({\displaystyle \genfrac{}{}{0pt}{}{ϵ,\frac{3}{2}2ϵ}{13ϵ}}|1c^2\right)`$
$`+e^{i\pi ϵ}{\displaystyle \frac{\left(1\right)_ϵ\left(1\right)_{2ϵ}}{\left(2\right)_{3ϵ}}}\left(1c^2\right)^ϵ`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,\frac{3}{2}2ϵ}{23ϵ}}\right|1c^2)\}.`$ (C.207)
The two hypergeometric functions are now combined into a single hypergeometric functions, which yields

$`{\displaystyle _Q}{\displaystyle \frac{2𝐩𝐫+r^2}{Q^2r^2\left(R_0^2+r^2c^2\right)}}_x`$ $`=`$ $`{\displaystyle \frac{1}{\left(4\pi \right)^2}}\left(e^\gamma \mu ^2\right)^ϵp^{2ϵ}{\displaystyle \frac{\left(1\right)_ϵ}{ϵ^2}}`$ (C.208)
$`\times \{{\displaystyle \frac{\left(\frac{3}{2}\right)_ϵ\left(1\right)_ϵ}{\left(\frac{3}{2}\right)_{2ϵ}}}+e^{i\pi ϵ}{\displaystyle \frac{\left(1\right)_ϵ\left(2\right)_{2ϵ}}{\left(2\right)_{3ϵ}}}(1c^2)^ϵ`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{ϵ,\frac{3}{2}2ϵ}{23ϵ}}\right|1c^2)\}.`$
Averaging over $`c`$, yields

$`{\displaystyle _Q}{\displaystyle \frac{2𝐩𝐫+r^2}{Q^2r^2\left(R_0^2+r^2c^2\right)}}_{c,x}={\displaystyle \frac{1}{\left(4\pi \right)^2}}\left(e^\gamma \mu ^2\right)^ϵp^{2ϵ}`$
$`\times {\displaystyle \frac{1}{ϵ^2}}{\displaystyle \frac{\left(1\right)_ϵ\left(1\right)_ϵ\left(\frac{3}{2}\right)_ϵ}{\left(\frac{3}{2}\right)_{2ϵ}}}\left[1+e^{i\pi ϵ}{\displaystyle \frac{\left(1\right)_{2ϵ}}{\left(1\right)_ϵ^2}}\right].`$ (C.209)
Expansion in powers of $`ϵ`$, yields Eq. (C.194).

### C.4 Hypergeometric functions

The generalized hypergeometric function of type $`{}_{p}{}^{}F_{q}^{}`$ is an analytic function of one variable with $`p+q`$ parameters. In our case, the parameters are functions of $`ϵ`$, so the list of parameters sometimes gets lengthy and the standard notation for these functions becomes cumbersome. We therefore introduce a more concise notation:

$$F\left(\genfrac{}{}{0pt}{}{\alpha _1,\alpha _2,\mathrm{},\alpha _p}{\beta _1,\mathrm{},\beta _q}|z\right){}_{p}{}^{}F_{q}^{}(\alpha _1,\alpha _2,\mathrm{},\alpha _p;\beta _1,\mathrm{},\beta _q;z).$$
(C.210)
The generalized hypergeometric function has a power series representation:

$$F\left(\genfrac{}{}{0pt}{}{\alpha _1,\alpha _2,\mathrm{},\alpha _p}{\beta _1,\mathrm{},\beta _q}|z\right)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{\left(\alpha _1\right)_n\left(\alpha _2\right)_n\mathrm{}\left(\alpha _p\right)_n}{\left(\beta _1\right)_n\mathrm{}\left(\beta _q\right)_nn!}z^n,$$
(C.211)
where $`\left(a\right)_b`$ is Pochhammer’s symbol:

$$\left(a\right)_b=\frac{\mathrm{\Gamma }\left(a+b\right)}{\mathrm{\Gamma }\left(a\right)}.$$
(C.212)
The power series converges for $`\left|z\right|<1`$. For $`z=1`$, it converges if $`\mathrm{Re}s>0`$, where

$$s=\underset{i=1}{\overset{p1}{}}\beta _i\underset{i=1}{\overset{p}{}}\alpha _i.$$
(C.213)
The hypergeometric function of type $`{}_{p+1}{}^{}F_{q+1}^{}`$ has an integral representation in terms of the hypergeometric function of type $`{}_{p}{}^{}F_{q}^{}`$:

$`{\displaystyle _0^1}𝑑tt^{\nu 1}\left(1t\right)^{\mu 1}F\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha _1,\alpha _2,\mathrm{},\alpha _p}{\beta _1,\mathrm{},\beta _q}}|tz\right)`$
$`={\displaystyle \frac{\mathrm{\Gamma }\left(\mu \right)\mathrm{\Gamma }\left(\nu \right)}{\mathrm{\Gamma }\left(\mu +\nu \right)}}F\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha _1,\alpha _2,\mathrm{},\alpha _p,\nu }{\beta _1,\mathrm{},\beta _q,\mu +\nu }}|z\right).`$ (C.214)
If a hypergeometric function has an upper and lower parameter that are equal, both parameters can be deleted:

$$F\left(\genfrac{}{}{0pt}{}{\alpha _1,\alpha _2,\mathrm{},\alpha _p,\nu }{\beta _1,\mathrm{},\beta _q,\nu }|z\right)=F\left(\genfrac{}{}{0pt}{}{\alpha _1,\alpha _2,\mathrm{},\alpha _p}{\beta _1,\mathrm{},\beta _q}|z\right).$$
(C.215)
The simplest hypergeometric function is the one of type $`{}_{1}{}^{}F_{0}^{}`$. It can be expressed in an analytic form:

$${}_{1}{}^{}F_{0}^{}(\alpha ;;z)=(1z)^\alpha .$$
(C.216)
The next simplest hypergeometric functions are those of type $`{}_{2}{}^{}F_{1}^{}`$. They satisfy transformation formulas that allow an $`{}_{2}{}^{}F_{1}^{}`$ with argument $`z`$ to be expressed in terms of an $`{}_{2}{}^{}F_{1}^{}`$ with argument $`z/\left(z1\right)`$ or as a sum of two $`{}_{2}{}^{}F_{1}^{}`$’s with arguments $`1z`$ or $`1/z`$ or $`1/\left(1z\right)`$. The hypergeometric functions of type $`{}_{2}{}^{}F_{1}^{}`$ with argument $`z=1`$ can be evaluated analytically in terms of gamma functions:

$$F\left(\genfrac{}{}{0pt}{}{\alpha _1,\alpha _2}{\beta _1}|1\right)=\frac{\mathrm{\Gamma }\left(\beta _1\right)\mathrm{\Gamma }\left(\beta _1\alpha _1\alpha _2\right)}{\mathrm{\Gamma }\left(\beta _1\alpha _1\right)\mathrm{\Gamma }\left(\beta _1\alpha _2\right)}.$$
(C.217)
The hypergeometric function of type $`{}_{3}{}^{}F_{2}^{}`$ with argument $`z=1`$ can be expressed as a $`{}_{3}{}^{}F_{2}^{}`$ with argument $`z=1`$ and different parameters 3F2 :

$`F\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha _1,\alpha _2,\alpha _3}{\beta _1,\beta _2}}|1\right)={\displaystyle \frac{\mathrm{\Gamma }\left(\beta _1\right)\mathrm{\Gamma }\left(\beta _2\right)\mathrm{\Gamma }\left(s\right)}{\mathrm{\Gamma }\left(\alpha _1+s\right)\mathrm{\Gamma }\left(\alpha _2+s\right)\mathrm{\Gamma }\left(\alpha _3\right)}}`$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{\beta _1\alpha _3,\beta _2\alpha _3,s}{\alpha _1+s,\alpha _2+s}}|1\right),`$ (C.218)
where $`s=\beta _1+\beta _2\alpha _1\alpha _2\alpha _3`$. If all the parameters of a $`{}_{3}{}^{}F_{2}^{}`$ are integers and half-odd-integers, this identity can be used to obtain equal numbers of half-odd-integers among the upper and lower parameters. If the parameters of a $`{}_{3}{}^{}F_{2}^{}`$ reduce to integers and half-odd-integers in the limit $`ϵ0`$ , the use of this identity simplifies the expansion of the hypergeometric functions in powers of $`ϵ`$ .

The most important integration formulas involving $`{}_{2}{}^{}F_{1}^{}`$ hypergeometric functions is (C.214) with $`p=2`$ and $`q=1`$. Another useful integration formula is

$`{\displaystyle _0^1}𝑑tt^{\nu 1}\left(1t\right)^{\mu 1}F\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha _1,\alpha _2}{\beta _1}}|{\displaystyle \frac{t}{1t}}z\right)`$
$`={\displaystyle \frac{\mathrm{\Gamma }\left(\mu \right)\mathrm{\Gamma }\left(\nu \right)}{\mathrm{\Gamma }\left(\mu +\nu \right)}}F\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha _1,\alpha _2,\nu }{\beta _1,1\mu }}|z\right)`$
$`+{\displaystyle \frac{\mathrm{\Gamma }\left(\alpha _1+\mu \right)\mathrm{\Gamma }\left(\alpha _2+\mu \right)\mathrm{\Gamma }\left(\beta _1\right)\mathrm{\Gamma }\left(\mu \right)}{\mathrm{\Gamma }\left(\alpha _1\right)\mathrm{\Gamma }\left(\alpha _2\right)\mathrm{\Gamma }\left(\beta _1+\mu \right)}}\left(z\right)^\mu `$
$`\times F\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha _1+\mu ,\alpha _2+\mu ,\nu +\mu }{\beta _1+\mu ,1+\mu }}|z\right).`$ (C.219)
This is derived by first inserting the integral representation for $`{}_{2}{}^{}F_{1}^{}`$ in (C.214) with integration variable $`t^{}`$ and then evaluating the integral over $`t`$ to get a $`{}_{2}{}^{}F_{1}^{}`$ with argument $`1+t^{}z`$. After using a transformation formula to change the argument to $`t^{}z`$, the remaining integrals over $`t^{}`$ are evaluated using (C.214) to get $`{}_{3}{}^{}F_{2}^{}`$’s with arguments $`z`$.

For the calculation of two-loop thermal integrals involving HTL averages, we require the expansion in powers of $`ϵ`$ for hypergeometric functions of type $`{}_{p}{}^{}F_{p1}^{}`$ with argument 1 and parameters that are linear in $`ϵ`$. If the power series representation (C.211) of the hypergeometric function is convergent at $`z=1`$ for $`ϵ=0`$, this can be accomplished simply by expanding the summand in powers of $`ϵ`$ and then evaluating the sums. If the power series is divergent, we must make subtractions on the sum before expanding in powers of $`ϵ`$. The convergence properties of the power series at $`z=1`$ is determined by the variable $`s`$ defined in (C.213). If $`s>0`$, the power series converges. If $`s0`$ in the limit $`ϵ0`$, only one subtraction is necessary to make the sum convergent:

$`F\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha _1,\alpha _2,\mathrm{},\alpha _p}{\beta _1,\mathrm{},\beta _{p1}}}|1\right)={\displaystyle \frac{\mathrm{\Gamma }\left(\beta _1\right)\mathrm{}\mathrm{\Gamma }\left(\beta _{p1}\right)}{\mathrm{\Gamma }\left(\alpha _1\right)\mathrm{\Gamma }\left(\alpha _2\right)\mathrm{}\mathrm{\Gamma }\left(\alpha _p\right)}}\zeta \left(s+1\right)`$
$`+{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}({\displaystyle \frac{\left(\alpha _1\right)_n\left(\alpha _2\right)_n\mathrm{}\left(\alpha _p\right)_n}{\left(\beta _1\right)_n\mathrm{}\left(\beta _q\right)_nn!}}`$
$`{\displaystyle \frac{\mathrm{\Gamma }\left(\beta _1\right)\mathrm{}\mathrm{\Gamma }\left(\beta _{p1}\right)}{\mathrm{\Gamma }\left(\alpha _1\right)\mathrm{\Gamma }\left(\alpha _2\right)\mathrm{}\mathrm{\Gamma }\left(\alpha _p\right)}}(n+1)^{s1}).`$ (C.220)
If $`s1`$ in the limit $`ϵ0`$, two subtractions are necessary to make the sum convergent:

$`F\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha _1,\alpha _2,\mathrm{},\alpha _p}{\beta _1,\mathrm{},\beta _{p1}}}|1\right)={\displaystyle \frac{\mathrm{\Gamma }\left(\beta _1\right)\mathrm{}\mathrm{\Gamma }\left(\beta _{p1}\right)}{\mathrm{\Gamma }\left(\alpha _1\right)\mathrm{\Gamma }\left(\alpha _2\right)\mathrm{}\mathrm{\Gamma }\left(\alpha _p\right)}}`$
$`\times [\zeta (s+1)+t\zeta (s+2)]+{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}({\displaystyle \frac{\left(\alpha _1\right)_n\left(\alpha _2\right)_n\mathrm{}\left(\alpha _p\right)_n}{\left(\beta _1\right)_n\mathrm{}\left(\beta _q\right)_nn!}}`$
$`{\displaystyle \frac{\mathrm{\Gamma }\left(\beta _1\right)\mathrm{}\mathrm{\Gamma }\left(\beta _{p1}\right)}{\mathrm{\Gamma }\left(\alpha _1\right)\mathrm{\Gamma }\left(\alpha _2\right)\mathrm{}\mathrm{\Gamma }\left(\alpha _p\right)}}[(n+1)^{s1}+t(n+1)^{s2}]),`$
where $`t`$ is given by

$$t=\underset{i=1}{\overset{p}{}}\frac{\left(\alpha _i1\right)\left(\alpha _i2\right)}{2}\underset{i=1}{\overset{p1}{}}\frac{\left(\beta _i1\right)\left(\beta _i2\right)}{2}.$$
(C.222)
The expansion of a $`{}_{p}{}^{}F_{p1}^{}`$ hypergeometric function in powers of $`ϵ`$ is particularly simple if in the limit $`ϵ0`$ all its parameters are integers or half-odd-integers, with equal numbers of half-odd-integers among the upper and lower parameters. If the power series representation for such a hypergeometric function is expanded in powers of $`ϵ`$, the terms in the summand will be rational functions of $`n`$, possibly multiplied by factors of the polylogarithm function $`\psi \left(n+a\right)`$ or its derivatives. The terms in the sums can often be simplified by using the obvious identity

$$\underset{n=0}{\overset{\mathrm{}}{}}\left[f\left(n\right)f\left(n+k\right)\right]=\underset{i=0}{\overset{k1}{}}f\left(i\right).$$
(C.223)
The sums over $`n`$ of rational functions of $`n`$ can be evaluated by applying the partial fraction decomposition and then using identities such as

$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{1}{n+a}}{\displaystyle \frac{1}{n+b}}\right)`$ $`=`$ $`\psi \left(b\right)\psi \left(a\right),`$ (C.224)
$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{\left(n+a\right)^2}}`$ $`=`$ $`\psi ^{}\left(a\right).`$ (C.225)
The sums of polygamma functions of $`n+1`$ or $`n+\frac{1}{2}`$ divided by $`n+1`$ or $`n+\frac{1}{2}`$ can be evaluated using

$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\psi \left(n+1\right)}{n+1}}{\displaystyle \frac{\mathrm{log}\left(n+1\right)}{n+1}}\right)={\displaystyle \frac{1}{2}}\gamma ^2{\displaystyle \frac{\pi ^2}{12}}\gamma _1,`$
(C.226)
$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\psi \left(n+1\right)}{n+\frac{1}{2}}}{\displaystyle \frac{\mathrm{log}\left(n+1\right)}{n+1}}\right)={\displaystyle \frac{1}{2}}\left(\gamma +2\mathrm{log}2\right)^2`$
$`+{\displaystyle \frac{\pi ^2}{12}}\gamma _1,`$ (C.227)
$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\psi \left(n+\frac{1}{2}\right)}{n+1}}{\displaystyle \frac{\mathrm{log}\left(n+1\right)}{n+1}}\right)={\displaystyle \frac{1}{2}}\gamma ^24\mathrm{log}2+2\mathrm{log}^22`$
$`{\displaystyle \frac{\pi ^2}{12}}\gamma _1,`$ (C.228)
$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{\psi \left(n+\frac{1}{2}\right)}{n+\frac{1}{2}}}{\displaystyle \frac{\mathrm{log}\left(n+1\right)}{n+1}}\right)={\displaystyle \frac{1}{2}}\left(\gamma +2\mathrm{log}2\right)^2`$
$`{\displaystyle \frac{\pi ^2}{4}}\gamma _1,`$ (C.229)
where $`\gamma _1`$ is Stieltje’s first gamma constant defined in (B.17). The sums of polygamma functions of $`n+1`$ or $`n+\frac{1}{2}`$ can be evaluated using

$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left(\psi \left(n+1\right)\mathrm{log}\left(n+1\right)+{\displaystyle \frac{1}{2\left(n+1\right)}}\right)`$
$`={\displaystyle \frac{1}{2}}+{\displaystyle \frac{1}{2}}\gamma {\displaystyle \frac{1}{2}}\mathrm{log}\left(2\pi \right),`$ (C.230)
$`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left(\psi \left(n+\frac{1}{2}\right)\mathrm{log}\left(n+1\right)+{\displaystyle \frac{1}{n+1}}\right)`$
$`={\displaystyle \frac{1}{2}}\gamma \mathrm{log}2{\displaystyle \frac{1}{2}}\mathrm{log}\left(2\pi \right).`$ (C.231)
We also need the expansions in $`ϵ`$ of some integrals of $`{}_{2}{}^{}F_{1}^{}`$ hypergeometric functions of $`y`$ that have a factor of $`\left|12y\right|`$. For example, the following 2 integrals are needed to obtain (C.166):

$`{\displaystyle _0^1}𝑑yy^{2ϵ}\left(1y\right)^{1+ϵ}\left|12y\right|F\left({\displaystyle \genfrac{}{}{0pt}{}{1ϵ,ϵ}{3ϵ}}|y\right)`$
$`={\displaystyle \frac{1}{6}}+\left({\displaystyle \frac{2}{9}}+{\displaystyle \frac{4}{9}}\mathrm{log}2\right)ϵ,`$ (C.232)
$`{\displaystyle _0^1}𝑑yy^{1+ϵ}\left(1y\right)^{1+ϵ}\left|12y\right|F\left({\displaystyle \genfrac{}{}{0pt}{}{2+2ϵ,1+ϵ}{2+3ϵ}}|y\right)`$
$`={\displaystyle \frac{1}{4}}+\left({\displaystyle \frac{7}{12}}+{\displaystyle \frac{2}{3}}\mathrm{log}2\right)ϵ.`$ (C.233)
These integrals can be evaluated by expressing them in the form

$`{\displaystyle _0^1}𝑑yy^{\nu 1}\left(1y\right)^{\mu 1}\left|12y\right|F\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha _1,\alpha _2}{\beta _1}}|y\right)`$
$`={\displaystyle _0^1}𝑑yy^{\nu 1}\left(1y\right)^{\mu 1}\left(2y1\right)F\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha _1,\alpha _2}{\beta _1}}|y\right)`$
$`+\mathrm{\hspace{0.33em}2}{\displaystyle _0^{\frac{1}{2}}}𝑑yy^{\nu 1}\left(1y\right)^{\mu 1}\left(12y\right)F\left({\displaystyle \genfrac{}{}{0pt}{}{\alpha _1,\alpha _2}{\beta _1}}|y\right).`$
The evaluation of the first integral on the right side gives $`{}_{3}{}^{}F_{2}^{}`$ hypergeometric functions with argument 1. The integrals from 0 to $`\frac{1}{2}`$ can be evaluated by expanding the power series representation (C.211) of the hypergeometric function in powers of $`ϵ`$. The resulting series can be summed analytically and then the integral over $`y`$ can be evaluated.