Cold quark matter at N3LO: soft contributions

High-order perturbative calculations for thermodynamic quantities in QCD are complicated by the physics of dynamical screening that affects the soft, long-wavelength modes of the system. Here, we provide details for the evaluation of this soft contribution to the next-to-next-to-next-to-leading order (N3LO) pressure of high-density, zero-temperature quark matter (QM), complementing our accompanying paper in arXiv:2103.05658. Our calculation requires the determination of the pressure of the hard-thermal-loop (HTL) effective theory to full two-loop order at zero temperature, which we go through in considerable detail. In addition to this, we comprehensively discuss the structure of the weak-coupling expansion of the QM pressure, and lay out a roadmap towards the evaluation of the contributions missing from a full N3LO result for this quantity.


I. INTRODUCTION
Determining the equation of state (EOS) of quantum-chromodynamic matter in extreme conditions using perturbation theory is a longstanding challenge almost as old as quantum chromodynamics (QCD) itself (see e.g. [2] for a review). In the case of high-temperature quark-gluon plasma (QGP), the calculation has reached a partial next-to-next-to-next-toleading order (N 3 LO) level [3][4][5]. At such high orders, a complication in perturbative calculations arises from the emergence of collective phenomena at long wavelengths, most importantly the physics of dynamical in-medium screening. To address this, all-loop-order resummations must be performed in order to reach a fixed order in the strong coupling constant α s .
At high temperatures T , reaching the partial N 3 LO accuracy was made possible on one hand by technical advances in the evaluation of multi-loop sum-integrals [6,7] and on the other hand by the seminal works of Kajantie et al. [3-5, 8, 9], where a resummation of soft screened modes of momentum scales α 1/2 s T and α s T was performed using the dimensionally reduced effective theories electrostatic QCD (EQCD) and magnetostatic QCD (MQCD) [10][11][12]. These calculations left only the contribution of the hard momentum scale πT missing from the full N 3 LO EOS of hot QGP, which constitutes a conceptually simple but technically very demanding challenge.
Screening phenomena closely analogous to those encountered at high temperatures appear also in the context of dense and cold quark matter (QM) [13,14], where phenomenological motivation stems from model-independent studies of the neutron-star matter EoS [15,16].
Here, the last fully completed order in perturbation theory dates back to the seminal papers of Freedman and McLerran [17,18], who determined the EOS to N 2 LO accuracy. At this level, the calculation becomes sensitive to the physics of screening, which these authors addressed through an all-loop-order diagrammatic resummation. The framework of dimensional reduction is unavailable at low temperatures, and challenges related to extending this resummation to higher orders have so far prevented bringing the EOS of cold QM to the same level of perturbative accuracy as its nonzero-T counterpart, although some progress in this direction has recently been achieved in [19,20]. In the present paper, complementing an accompanying letter [1], we finally perform this resummation using the Hard-Thermal-Loop (HTL) effective theory, determining the soft contributions to the EOS up to and including the N 3 LO order. While Ref.
[1] concentrates on an in-depth analysis of the result, here we provide extensive details of the technical aspects of the calculation, and in addition discuss the computations needed to determine the last contributions missing from a full N 3 LO result for the EOS of cold QM.
The physical picture behind perturbative calculations at high densities is as follows. In a medium characterized by a large quark chemical potential µ q and zero temperature, cold QM contains a filled Fermi sea of quarks from zero momentum up to the scale µ q . 1 The free Fermi pressure of this system of quarks forms the leading-order (LO) description of the pressure p of cold QM and scales as µ 4 q in the case of massless quarks. 2 While there are no on-shell gluons in the medium, off shell gluons are present because the quarks are color charged. Interactions between the quarks and gluons in QCD lead to corrections to this LO pressure as a function of the strong coupling constant α s .
Because of this Fermi sea of quarks, the propagation of both quarks and gluons through cold QM becomes modified. Low-momentum quarks are Pauli blocked and cannot propagate, as those states are filled by the medium. Thus, these low-momentum quarks do not contribute to higher-order loop corrections to p, leaving the scale µ q (dubbed "hard") as the only relevant scale for quarks. For gluons, the picture is more complicated and involves two different approximations that can be used in different regions of momentum space: the naive loop expansion and the HTL expansion (see Fig. 1). Hard, short-wavelength modes can be treated similarly to the quarks in a naive loop expansion, while the soft, long-wavelength gluons become qualitatively modified by the medium and require resummations of arbitrary numbers of one-quark-loop insertions within calculations. These modifications lead to, e.g., nonanalyticities ln α s in the weak-coupling expansion in the pressure of cold QM.
The rest of this Introduction is organized as follows. In Sec. I A, we introduce the naive loop expansion and the HTL expansion and motivate their respective regions of validity.
After this, in Sec. I B we explain how to power-count the contributions of the resummed soft gluons. In Sec. I C, we then discuss the analytic structure of the different contributions to the pressure of cold QM, proceeding from LO to N 3 LO. Finally, in Sec. I D we explain what precisely is computed in the present article, and walk the reader through the overall structure 1 In this section, we shall describe the situation for a single quark flavor, for simplicity. 2 Note that at high density, quarks on the Fermi surface undergo pairing through attractive channels of gluonic interactions, leading to a different ground state [21]. However, these effects do not enter at any finite order in the weak-coupling expansion. At sufficiently high densities, the pairing gap ∆, which only depends on µ in a mild way, becomes small in comparison to the chemical potentials, and the pairing contributions to the pressure become suppressed by a factor ∆ 2 /µ 2 . of the paper.

A. Two expansions for gluons in cold quark matter
Whether gluonic propagation is qualitatively modified by scattering from hard quarks in loop corrections depends on the magnitude of the propagating gluonic momentum [22]. This can be seen most clearly from the dispersion relation of the gluonic modes with momenta P , which is schematically of the form P 2 + Π(P ) = 0. (1) Here we will work consistently in a Euclidean framework, and Π(P ) is a generic component of the Euclidean gluon polarization tensor. This tensor is parametrically of the order of the square of the in-medium effective mass scale m E , related to the one-loop Debye mass. For a single massless quark in d = 3 spatial dimensions it has the value m 2 E = (2/π)α s µ 2 q . 3 If the free part of Eq. (1) parametrically exceeds the interaction part, the interactions can be treated as perturbations to the free propagation of gluons and be dealt with using a naive perturbative (loop) expansion. We can see that this occurs for |P | m E .
If, on the other hand, the gluon has momentum |P | m E (dubbed "soft"), then its propagation is qualitatively modified. In particular, generic low-momentum gluonic excitations require a nonzero excitation energy proportional to m E . 4 This behavior arises because the gluonic self energy Π µν (P ) has a nonzero |P | → 0 limit lim |P |→0 Π µν (P ) ≡ Π µν HTL (P ) = 0 (2) 3 This effective mass scale is related to the asymptotic HTL mass [23] m ∞ by m 2 E = 2m 2 ∞ . In the case of multiple quark flavors f , this effective mass scale becomes m 2 E = (2α s /π) f µ 2 f . 4 Note that for some special directions of the Euclidean four-vector P , these excitations may still be massless. At high temperature, unscreened magnetic gluons lead to the generation of a further ultrasoft mass scale, but this is not the case in cold QM, as soft gluons are not Bose enhanced (see below). where we have suppressed the color indices and defined a unit four-vectorP = P/|P | in the direction of P , a notation we shall use prominently in this work. Here, we have also identified the HTL self energy, which is the low-momentum limit of the full self energy. It is important to note that in cold QM, only the quark loops contribute to this HTL self energy: hard gluon and ghost loops are not directly populated by the medium (i.e., there are no on-shell gluons present) and thus do not differ from their vacuum values (i.e. depend on µ q ) at LO. It then follows that only the single quark loops must be resummed into the gluonic propagators, as all other corrections can be naively expanded around that limit. This resummation is depicted in Fig. 2.
The full kinematics of the one-loop quark contribution to the gluon polarization tensor results in a Π µν (P ) that depends not only onP , but also on the magnitude of the gluonic momentum P . However, if this magnitude is parametrically less than the hard scale, |P | µ q , one may systematically expand the polarization tensor in powers of the ratio of the external gluon momentum and the internal quark loop momentum, |P |/|Q|, which constitutes part of the HTL framework as described by Braaten and Pisarski [24]. In addition to this modified propagation, the interaction between soft gluons also becomes modified in cold QM, as n soft gluons interacting through a quark loop with momentum |Q| µ q enters at the same order as the bare coupling between n soft gluons (if it exists). For example, for soft gluons the interaction shown in the right panel of Fig. 2 enters at the same order as the bare three-gluon coupling. However, both quark and ghost fields remain unresummed in Euclidean spacethe quark propagators are protected in the infrared by the nonzero chemical potential, while the ghosts are known not to develop a thermal mass [25].
In summary, there are two different approximations that can be made for gluons in cold QM, depending on the magnitude of their momenta P : if |P | µ q , the HTL expansion becomes valid, and if |P | m E , the naive loop expansion becomes valid. This is shown pictorially in Fig. 1. As demonstrated in Ref. [20], the integration region m E |P | µ q , where both approximations are valid (dubbed "semisoft" in Ref. [20]) leads to a logarithm of the coupling: Here, the notation · P indicates an average over all four-dimensional Euclidean angles. Since such a logarithm can arise from any integral over a resummed gluonic momentum, at higher orders, where multiple resummed gluons may contribute, we should expect to find contributions to the pressure containing factors of ln n α s , where n is the number of resummed gluonic momenta in a given resummed diagram. At N 3 LO, this leading logarithm was computed already in Ref. [20].
We now address the question of how one should power count such resummed, soft gluons, to see where they first contribute.

B. Power counting the soft contributions
As per the discussion in Ref. [20], the soft gluons occuring in loop corrections and that require resummation are phase-space supressed by the integration measure where the two effective vertices are represented by Γ(K, P, R) 2 in the numerator. For soft momenta, these vertices scale in the same way as the bare vertices, leading to Γ(K, P, R) 2 ∼ m 2 E . Similarly, for soft P, K, R, the denominator scales as m 4 E , leading to the final power counting shown above. We thus see that this is a perturbative correction to the self energy Π(K) ∼ m 2 E . The case of adding a gluon through a four-point vertex is identical, as for the purposes of power counting a single 4g HTL vertex behaves similarly to two 3g HTL vertices. We therefore see that these soft corrections are indeed perturbative. A full account of the form of these effective vertices, and our notation for the double lines and blob vertices for the HTL diagrams are introduced in detail in Appendix B.
The fact that interactions between soft gluons are perturbative means that we can systematically improve the soft sector using a loop expansion within the HTL theory. This situation is qualitatively different from that encountered in high-temperature QGP, where the presence of "ultrasoft" gluonic momenta of order α s T famously leads to the Linde problem and the emergence of fundamentally nonperturbative contributions to the pressure at O(α 3 s ) (see e.g. Ref. [2] for a discussion of this subtle topic).
Let us now make one brief remark about regularization. The cutoff description implicitly used so far in these discussions (and in Ref. [20]) is very convenient for identifying the physical sources of the logarithms. However, in detailed computations (especially at higher orders) such an approach has drawbacks. In particular, once there are multiple soft gluons, performing the entire computation with a cutoff is very cumbersome. Thus, at this point we make the choice to use dimensional regularization to regulate not only the hard UV divergences arising in the full theory but also the intermediate ones arising from different kinematic regions (which we will introduce momentarily), once those regions are separated.
This will provide a much more streamlined framework for self-consistently determining all the different contributions to the pressure at higher orders.
C. Computing the pressure of cold quark matter Following the logic from above, we deduce the following structure for the pressure of cold QM, valid up to and including the N 3 LO terms: Here, p FD is the pressure of a free Fermi gas of quarks, while the other terms arise from interaction corrections among or across modes of different types. Terms on first line arise from hard modes and can be computed through a naive loop expansion in full QCD. Terms on the second line arise from soft modes and their interactions, and can be determined within the HTL theory. Finally, the remaining term on the third line arises from interactions between the soft and hard modes and requires a partial HTL resummation.
Due to the ambiguous semisoft momentum range m E P µ q , the splitting between the different kinematic regions is not unique. This ambiguity leads to ultraviolet (UV) divergences within the p s i that cancel against corresponding infrared (IR) divergences within the p h i (and mixed UV-IR divergences in p m 3 at N 3 LO). This cancellation will be further remarked on briefly below. The ambiguity also makes these coefficients dependent on a factorization scale Λ h , which arises from the dimensionally regularized integration measure in our case, and which will be canceled when summing over the different kinematic contributions at a given order. In terms of the factorization scale, the divergences also lead to expressions of the form ln(Λ h /m E ) and ln(µ q /Λ h ) from the UV limit of the soft sector and the IR limit of the hard sector, respectively. As the Λ h dependence cancels in the sum over all kinematic contributions, these logarithms will generate precisely the ln(µ q /m E ) ∼ ln α s terms discussed in Sec. I A. It is important to note here that the factorization scale Λ h is not a momentum cut-off between the soft and hard sectors of the theory, and thus need not lie between the scales m E and µ q . This becomes relevant when analyzing the behavior of the result, and is further discussed in Ref. [1].
To discuss the structure of the terms in Eq. (6) in more detail, we find it useful to further classify the diagrams contributing to the different coefficients p n i . This further classification is based on the HTL limits of the different diagrams which are obtained by taking soft kinematics of the all the gluon lines. Diagrammatically this is reached by (i) resumming all gluon lines and (ii) contracting all quark loops into points and absorbing them into propagators and HTL vertices (see a sample illustration of this process of "HTL resummation" in Fig. 3). Concretely, we classify the contributions by the loop order of the resulting HTL diagram after taking fully soft kinematics. 5 In general, for 0 < j < i, an (i + 1)-loop hard diagram is classified as a part of p h i,j if its fully soft limit is a j-loop HTL diagram. The remaining (i + 1)-loop hard diagrams, which lead to i-loop HTL diagrams, are classified as This leads to the decomposition of the N 2 LO and N 3 LO terms and where the second subscript denotes the number of hard gluon (or ghost) loop momenta, except for the above-mentioned exceptional p h i,0 terms. We show the relation between these contributions and the fully resummed HTL diagrams in Fig. 4. Organizing the terms in this fashion guarantees that the sum of all contributions in a given column is independent of the factorization scale as well as associated divergences. 6 Next, we discuss in detail each of the terms appearing in Eqs. (6)- (8). • LO: This is simply the free Fermi pressure, arising from a single diagram with full quark kinematics; • NLO: In this contribution, only a hard gluon contributes, since a soft gluon is phasespace suppressed. The quark loop requires the full kinematics.
At N 2 LO, there are contributions from both hard and soft gluons, corresponding to the coefficients p h 2 and p s 2 in Eq. (6), respectively. As the diagrams in p s 2 are fully resummed, the diagrams are IR safe and there is only one subclass of diagrams contributing to On the other hand, p h 2 can be further subdivided into two subclasses and grouped by their IR properties. Here, the hard p h 2,1 and soft p s 2,0 become by construction identical for semisoft gluon kinematics, and the diagram in p s 2,0 can in fact be generated from that in p h 2,1 by HTL resumming the gluon loop in the hard diagram. One might worry that this leads to a double counting of contributions, but this is not the case since at T = 0 the semisoft region of p h 2,1 or p s 2,0 gives rise to scalefree integrals that vanish in dimensional regularization: in the first case, the quark loops can be replaced with HTL self-energy insertions, in the second, the resummed gluon propagator is re-expanded; in both cases, the integration over the gluonic loop momentum is then scalefree in the semisoft region. 7 The p h 2,0 contribution is on the other hand IR safe and diagrammatically distinct. This is precisely the structure illustrated in the left panel of Fig. 4.
The N 2 LO contributions to the pressure of cold QM were first determined by Freedman and McLerran in Refs. [17,18] without the use of the HTL theory. In the modern language of dimensional regularization, the logarithmic contributions to the pressure arise solely from the subclasses p s 2,0 and p h 2,1 above, through the (1/ε) × ε terms therein, with ε ≡ (4 − D)/2 and D standing for the spacetime dimensionality. Additionally, one finds 1/ε terms cancelling between these p s 2,0 and p h 2,1 subclasses. In the case of p s 2,0 , such a term arises from a UV divergence, while in p h 2,1 , it arises from an IR divergence.

Classification of diagrams at N 3 LO
The organization of the N 3 LO contributions was displayed already in Eq. (8) and the right panel of Fig. 4 and is further visually summarized in Fig. 5.
Similar to the N 2 LO calculation, at N 3 LO there are multiple classes of contributions, i.e. those which arise from either two, one, or zero soft gluons, corresponding to the coefficients p s 3 , p m 3 , and p h 3 in Eq. (6), respectively. In this section, we give a more detailed account of all the subclasses making up these contributions, proceeding according to the HTL diagrams they are related to (i.e., the colored columns in Fig. 4).
The soft contribution is again fully resummed and IR safe and therefore forms only one subclass, namely These diagrams are intimately related to those in the mixed contribution 8 α 3 s p m 3,1 = The diagrams in p m 3,1 become identical to those in p s 3,0 in the semisoft region, while the diagrams of p s 3,0 can be generated from those in p m 3,1 by HTL resumming the one unresummed 8 Here and in the following diagrams, we assume an implied summation over the fermion and ghost directions in each loop, to reduce the number of diagrams shown. One can find the full list of contributions, with the correct symmetry factors, in Ref. [27].
gluon line. Furthermore, the kinematics of the resummed gluon line in this mixed contribution is soft, and hence one may expand in the small gluonic momentum. In this sense, these hard corrections within the mixed p m 3,1 contributions can be thought of as corrections to the HTL self energy.
Similarly, the diagrams in are intimately related to those in p m 3,1 above. Upon HTL resumming one of the gluon lines, one obtains the diagrams of p m 3,1 . Now we turn to the second column in the right panel of Fig. 4 related to the one-loop HTL diagrams. The one-loop diagram contributes at N 2 LO (namely, in p s 2,0 ), so to contribute at N 3 LO it has to be dressed with a hard quark line. 9 This comes about naturally when one of the momenta running in the HTL self energies becomes hard, giving rise to the diagram where we can again expand in the soft gluon momentum. However, since the leading-order term in the small-momentum expansion gives back the lower-order one-loop HTL result, we must instead use the NLO soft kinematics [28][29][30] to obtain the N 3 LO contribution.
The diagram in p m 3,0 is on the other hand related to the graph 10 For semisoft kinematics, the diagram in p m 3,0 becomes identical to the one in p h 3,1 . Finally, there remains one further subclass where no resummations are necessary, namely p h 3,0 which contains the remaining IR-safe four-loop diagrams containing a single quark loop, here with the full quark kinematics. It reads The hard gluon or ghost contributions already appear in p m 3,1 , and are related instead to the two-loop HTL diagrams. 10 It is worth noting that this diagram represents the leading large-N f behavior of the pressure, which has been determined in Ref. [31] in the high-T limit.
The set of all four-loop vacuum diagrams in QCD has previously been written down using a different organization in e.g. Ref. [27]. Comparing the two, we see that each diagram is correctly reproduced in our organization, and we have checked agreement with the symmetry factors as well.

D. What we compute in this work
In this paper, we determine the contribution p s 3 to the cold-QM pressure, defined in Eq. (6), which is equal to the fully soft subclass p s 3,0 at N 3 LO. While it does not constitute a complete new order in the weak-coupling expansion of the pressure, it amounts to a complete kinematic contribution that has furthermore been speculated to play a crucial role in the slow convergence of the quantity [32]. In addition, as we will show below, we can recover the known O(α 3 s ln 2 α s ) contribution to the pressure from this region alone.
In the computation, we use dimensional regularization and work in the limit of vanishing quark masses, which amounts to evaluating the full two-loop HTL pressure at zero temperature, without expanding in the in-medium effective mass scale. We note that there has been previous research on higher-order HTL thermodynamics [32][33][34][35][36]. However, these works all expand the HTL diagrams in powers of the effective mass, and so do not perform the full resummation that we need.
The general structure of the paper is as follows. In Sec. II, we introduce the setup, conventions, and machinery used in the calculation of p s 3 . We explain the power counting of all the contributions in detail and present notations that allow us to easily extract the UV-sensitive integral contributions. In Sec. III, we then explain our steps for evaluating these integral expressions, and display results for the different contributions along with many details, especially for the UV-sensitive terms. Finally, in Sec. IV, we present our final result for the pressure in the soft region.
We also discuss cross-checks of our computation, remark upon the sizes of different contributions, and provide a small outlook for the remainder of the full N 3 LO pressure. Following the main text is a large collection of appendices summarizing the Euclidean-space HTL framework used throughout this work, as well as the additional machinery that we have developed to tackle the computation. We have collected all of this into the appendices to aid future researchers who wish to use the Euclidean-space HTL framework in their work.

A. Starting expression and convention
The expression corresponding to the fully soft contribution to the cold-QM pressure that we will evaluate in this paper is [cf. Eq. (6)] with I 3g , I 4g , I gh labeling the diagrams repeated in Fig. 6. Here, g = √ 4πα s is the QCD gauge coupling, N c is the number of quark colors, and d A = N 2 c − 1 is the dimensionality of the adjoint representation of the SU(N c ) gauge group, or the number of gluons. We will perform this computation in dimensional regularization in D = d + 1 spacetime dimensions.
Eq. (20) is the expression for the two-loop HTL pressure examined also in Ref. [37], but in the zero temperature limit and at nonzero density. To evaluate this expression, we find it easier manipulate it using techniques that will be discussed below, rather than using the expressions in Ref. [37] as a starting point.
First, we write down explicit expressions for the resummed two-loop graphs under study.
Using the Feynman rules of Appendix A, we readily obtain These expressions form the starting point of our diagrammatic analysis. We have here introduced the following notations, which are discussed in more detail in the appendices: Metric and vector conventions: We work in Euclidean space with metric δ µν . We write the components of our four-vectors as K µ = (K 0 , k i ), where K 0 = K 0 , k i = k i for i = 1, 2, 3.
The scalar product between two Euclidean four-vectors K and P is given by K · P = K 0 P 0 + k · p, and we use the notation |k| 2 = k · k for the magnitudes of the spatial part of the momentum. We will also repeatedly use the notation for the unit vectors in the direction of K and k, respectively, the former already defined around Eq. (2).
Integration measures: The integration measures are defined in D = d + 1 space-time dimensions. The symmetric integral KP R we use frequently is defined as where Λ h is a factorization scale and the factor (e γ E /4π) (4−D) , with γ E the Euler-Mascheroni constant, is introduced so that one absorbs the ultraviolet (UV)-divergent part with a universal constant. The integral KP in Eq. (22) follows from KP R upon doing the trivial integration over R.
Propagator: The HTL-resummed gluon propagator D µν is defined in the covariant gauge as where the parameter ξ fixes the gauge and In this equation, and from this point on, we drop the HTL label on the HTL self energies, for brevity. It turns out that the computation is most efficient to perform in the ξ = 1 gauge, which we shall use throughout the rest of the text. In the class of covariant R ξ -gauges, the gauge parameter ξ only appears in the propagator, and we have explicitly checked that the expression in Eq. (20) is gauge-independent in the sense of being independent of the ξ-parameter, with the conclusion being supported by Ref. [37].
The standard projectors P µν T (K) and P µν L (K) used in Eq. (26) are defined in Sec. B, and the reader is advised to consult the appendix when necessary. We do repeat here the definition of the symmetric and transverse HTL self-energy tensor where the coefficient functions are given by The trace of the one-loop HTL self energy is defined to be m 2 E , so that m 2 E ≡ Π µµ (K). An alternative explicit definition of Π µν (K) is given in the aforementioned appendix.
Vertices: The effective three-(3g) and four-gluon (4g) vertices are obtained by adding the HTL loop to the bare vertex. We write these quantities as where the subscript "0" is always understood as referring to the bare quantity. The HTL vertices δΓ are only defined when the sum of all of their arguments is zero, and they have the property that they are totally symmetric in their indices and traceless in any pair of indices. Furthermore, the 3g HTL vertex is even (odd) under even (odd) permutations of the momenta P , Q, and R, while the 4g HTL vertex is even under all permutations of the momenta P , Q, R, and S. We note here that like the HTL self energies, the HTL vertices are proportional to m 2 E , and satisfy the generalized Ward identities We also note that the HTL vertices are independent of the gauge parameter ξ, even without fixing the gauge. The explicit integral expressions for the d-dimensional HTL vertices are given in Appendix B 2.

B. Isolating the UV-sensitive terms
Due to the two-loop structure of the quantity under consideration and the fact that the HTL self energies are independent of the magnitude of the gluonic momenta, we can deduce that the general structure of our final result for α 3 s p s 3 will be When this is expanded out fully in ε, our final result for the soft contributions to the pressure of cold QM becomes and contains both double and single-logarithmic terms of the ratio m E /Λ h . Note that, since we have resummed the soft sector, the expression is IR safe. Thus, the terms in Eq. (32) that enter with negative powers of ε arise from the UV. We are thus led to the following conclusion: in order to isolate the p −2 and p −1 parts of the pressure, we must isolate the UV behavior of the integrals Eqs. (21)- (23). To do this, we will first power count to isolate the power-law divergences in the UV (which do not contribute at all in dimensional regularization). This in turn will lead us to introduce a notation to isolate the terms identified by the power counting.

Power counting
In the UV, the bare and HTL parts of the vertices and propagators no longer enter at the same order in m E . Therefore, by expanding the vertices and unfolding the propagators, we will be able to isolate the UV-sensitive terms. Let us first begin by expanding out the vertices into their bare (0) and HTL (H) parts, and splitting I 3g and I 4g into the following pieces: with and where we have used a compact notation for the power of Π µν (K), and have used the simple form of the bare gluonic propagator in ξ = 1 gauge in the leading term, namely Taking a momentum K of the order of some large Λ m E , the expansion in Eq. (37) goes as which shows that we have an infinite number of contributions, with the higher terms becoming less important. On the other hand, when the momentum flowing through the vertex becomes soft K ∼ m E , the propagator simply scales as so that no expansion is possible.
Let us now use these scalings to power count one piece of Eq. (35). When K, P, R ∼ Λ, the term with two bare vertices scales as which shows that there is also a power-law divergence here, even though one of the propagators is still resummed. Observe that letting one of the lines become soft shifts the leading term in the expansion of the integral to a higher power in m E . Note also that in this region, there are subleading contributions from the vertices, in addition to the propagators. Finally, there is the soft-soft region. However, this region does not probe the UV, and so will not lead to any divergences at all.
Inspecting the above equations, we see that these power-law divergent terms arise from the first few leading terms in the expansion of the propagators, which motivates us to introduce a new notation in the following section. 11 Note that we are power counting in D = 4.

Peeling away the bare propagators
We have called the leading term in Eq. (37) (the bare propagator) ∆ µν 0 (K); let us extend this notation to label the other terms in the expansion as well: Here we make the identification [Π(K) 0 ] µν = δ µν to match the leading term. This allows us to write the expansion in Eq. (37) as We can now introduce the following notation for the resummed propagator with the n leading terms removed: Consequently, the D µν n (K) are still resummed expressions, while the ∆ µν n (K) are not. Note that both ∆ µν n (K) and D µν n (K) are D-dimensionally transverse for every n, and that the following relations hold for any n ≥ 1: Notice also the full propagator at the end of both lines of Eq. (47), and the fact that Eq. (48) is not a partial sum of an infinite series, but is exact.

Applying the new notation
Using this new notation, we can make the schematic expressions in Eqs. (42)-(43) more explicit. To this end, we rewrite I (0,0) 3g into the following form: We stress that this equation is not an expansion, but holds exactly. 12 Observe that the two terms on the second line of the above equation contain only ∆s, which have no mass scale. Therefore, these two terms are power-law divergent and thus vanish in dimensional regularization. These terms correspond precisely to the two power-law divergent terms in Eq. (42). The remaining power-law divergence in Eq. (43) is located in one particular part of the term containing D 2 on the third line of the above equation, as can be seen by counting the bare propagators. This divergence will not become explicit until we insert the explicit expression for the vertices and perform the contractions.
In addition to identifying the terms that contain only ∆s and are trivially power-law divergent in the UV, this notation also allows us to identify precisely the terms which are logarithmically UV-sensitive. We can do this as follows. From the scaling of the ∆ n and D n in Eqs. (50)-(51), we see that we can determine the scaling of a product of ∆ n s and D n s by simply summing the subscripts: if the sum is N then the product goes as m 2N E Λ −2(N +1) in the UV. However, we know that the logarithmically UV-sensitive terms (ones that lead to the O(ε −2 ) and O(ε −1 ) terms in the expansion) must scale as m 4 E times a dimensionless integral in the UV, as such a structure allows it to contribute through the whole UV tail of the integral. This gives us a powerful prescription for identifying the logarithmally UV-sensitive terms: we can simply sum the subscripts on the propagator terms (and include one m 2 E per HTL vertex if necessary) and see whether the result is m 4 E . We now carry out the above procedure for all of the terms in Eq. (34) as well as for I gh to peel away all the power-law UV divergences and to identify the logarithmically UV-sensitive 12 Note that, from now on, we will use this more compact expression for the vertices.
terms. We then find the following: for the parts of I 3g , for the parts of I 4g , and for I gh , Here, we have introduced the notation [I] ijk to mean the integral I with the propagators (In the case of fewer subscripts, the notation is the same, but replacing only those propagators appearing in I. ) We have additionally added a superscript "UV" to those terms which are logarithmically UV sensitive, and hence contribute to the O(ε −2 ) or O(ε −1 ) terms of the final answer. From this point on, we shall simply refer to them as the "UV terms". Consequently, those terms which do not have the UV labeling are UV finite and thus can be computed in d = 3 spatial dimensions. Most of these non-UV-sensitive terms cannot be analytically simplified much more, and so for the rest of this organizational section we shall not manipulate them. Their contributions are listed in Sec. III C.

Performing the contractions
In most of the UV terms, the Lorentz contractions are relatively straightforward and lead to some simplifications. The general procedure for evaluating all of them, except those in UV , which we will separately consider in Sec. III B, is as follows: 1. Substitute in the form of the bare propagator ∆ 0 in Eq. (39) to the expressions, which contracts some of the indices of the vertices together. 3. Use the symmetry of the integration measure KP R to permute K, P, R as necessary to further simplify the expressions.
Let us quickly illustrate this procedure in two cases, one with only bare vertices I (0,0) 3g UV 002 , and one with a 3g HTL vertex I (0,H) 3g UV 001 . In the first case, we have where from the first to the second line, we used step (1) above; from the second to the third line, we used step (2); and from the third to the fourth line, we used step (3) of our procedure. Finally, in the last line, we separated out a UV power divergence, which vanishes in dimensional regularization. This is precisely the power-divergent term identified in Eq. (43). Here, we have introduced a Lorentz-index-free notation, where we represent contractions as dots or traces depending on the structure: We will use this notation extensively through the rest of the work.
In the second case, we have where we used the same steps in the same order. Here, the generalized Ward identity was applied in the contraction of the bare 3g vertex and the 3g HTL vertex. Additionally, we also made use of the symmetries of the 3g HTL vertex discussed in Appendix B 2.
Now let us further massage the above expression by using (i) the fact that for any n (this follows from Eq. (47)), and (ii) the fact that, since P 2 ∆ µν 0 (P ) = δ µν , where in the last line we used the recursion relation in Eq. (62) immediately above. Thus, plugging in the above and rearranging, we find Note that the last term in this expression goes as m 6 E in the UV, and so can only contribute to the constant term.
Our final results for all the UV terms after the above procedure are as follows: for the parts of I 3g , for the parts of I 4g , and for I gh , The sum of the above terms is quite compact, as there are a few cancellations. The total reads: Here, we have factored out all magnitudes from the dot products, using our unit-vector notation introduced in Eq. (24). Note the very important point that the HTL self energy Π µν (K) inside of the propagators depends only onK, and not on the magnitude |K|. In light of this, we will from now on write Π µν (K) to make this clear.
At this point, we are ready to tackle performing the integrations of all the terms isolated and scalarized above.

III. THE COMPUTATION
After the manipulations performed in the previous section, evaluating p s 3 has been reduced to computing [I tot ] UV in Eq. (71) and the sum of the non-UV terms in Eqs. (53)-(58). Our general procedure for evaluating these integrals can be summarized as follows: 1. Rescale all the magnitudes of momentum variables by |K| → m E |K| etc. to pull out the m E dependence and make the integrands dimensionless.
2. Perform the trivial R integral to set R → −(K + P ), and change variables in the remaining K and P integrals to write them as integrals over the magnitude of the four-vectors |K| and |P |, and the remaining angles.
3. Further transform from the magnitudes of the momenta (|K|, |P |) to Euclidean polar (X, χ) coordinates, given by We introduce a shorthand notation s χ , c χ , as these particular functions will appear many times in our computation. In all cases, the radial X integral can be performed analytically (in general d, if necessary).
Let us now discuss the details of these operations.
Details of step (1): As part of this step, we will define dimensionless versions of all of the functions, denoted with tildes. For instance: where D n (K) and G X (K) contain Π µν (K) in place of Π µν (K). This rescaling also introduces a change in the integration measures; for example where we have defined a new integral (again, KP can be defined by performing the trivial integral over R) by adding m 2 E to the denominator of the factor out front. Thus, by factoring out m 8 E , and by simply "placing tildes on everything else", we can implement the desired rescaling of variables. In what follows, we shall assume that we have already done this procedure everywhere, and hence drop the explicit tildes for convenience. (2): After performing the trivial integral over R, which sets R → −(K + P ), we wish to change variables to radial and angular parts. To this end, we introduce the following notations. We first define the Euclidean spherical variables (|K|, Φ K ), where Φ K ∈ [0, π] is the D-dimensional polar angle defined by

Details of step
The dot product between the Euclidean four-vectors K and P can be written as where the magnitudes |K| and |P | factorize from the angular part w ≡ w(Φ K , Φ P , θ) with cos θ ≡k ·p. Note that this equation implies thatK ·P = w, which will be used extensively in the computation below. Now, in these new coordinates, the integration measure KP can be written where the angular part of the integration measure is and we have introduced a compact notation for the dimensionless prefactor The integrals over the azimuthal angles have been performed and included in C(d), as all of our integrands are independent of those angles. For reference, we note here that C(3) = 2/(2π) 6 , and that in d = 3, the integral Ω = π 2 /2.

Details of step (3):
Here, we stress that by moving to the (X, χ) coordinates, we only change the magnitudes |K|, |P |, and soK,P remain unchanged. This change of variables gives rise to a nontrivial Jacobian, which leads to Finally, note that since we only integrate over the first quadrant in (X, χ) space, we have s χ , c χ ≥ 0 everywhere in our integrals.
We are now in a position to evaluate our terms. We begin by identifying the independent integrals in [I tot ] UV : which combine as follows Note that we have kept the (d − 1) inside I D , since we want to be able to do the ε expansion for each of these integrals separately. Additionally, I E is not UV sensitive, as can be seen by counting the propagator subscripts as before. Therefore, we will now proceed as follows.
First, as the forms of I A -I D are so similar, and since they are all UV sensitive, we shall treat all of them together in Sec. III A below, and define I ABCD to be the corresponding part of Eq. (89): so that Second, as the structure in the remaining UV-sensitive term I F is so different from the others (in particular, the 3g HTL vertices cannot be removed completely by using the generalized Ward identities), we will evaluate it separately in Sec. III B. Finally, as all the remaining finite terms, including I E , can be numerically evaluated in d = 3, we treat them all in Sec. III C.
A. Integrals I A -I D After performing the steps (1)-(3) described in the previous section, we are left with integrals over the angles in Ω , the separate angle χ, and the radial coordinate X. As alluded to above, the X integral can be performed analytically in general d. In I A -I D only the following two master integrals in X appear: Here, I, J ∈ {T, L} label the polarization components of the HTL self energy. We will now proceed through the integrals one by one, mostly showing details only for I A .
1. Doing the integrals in general d To unpack the notation in I A , we use the definitions for the propagators D n and D in Eqs. (47) and (26), respectively. This leads to the following form for I A in the (X, χ) coordinates where C(d) was defined in Eq. (81). Using the first master integral Eq. (92), this becomes The integral over χ here can be performed analytically, but it is unwieldy, and so we choose the following approach instead. This integral contains a divergence in d = 3, from the region near χ = 0. Near χ = 0, the χ integrand is approximately This means that we can isolate the O(ε −2 ) part of I A by evaluating where in going to the second line, we used Tr [Π] = 1, which is true because we have factored out the mass scale m E , and have dropped the tildes on all expressions. We have also introduced the notation We will return to this χ ≈ 0 term in a bit.
The remaining integral over χ, which we define to be in the following, will then be finite: Note that this χ integral vanishes in d = 3, as the expression in braces equals 0. Thus, we only need the O(ε) piece of the χ integral, and the O(1/ε) piece of the divergent coefficient The full expression up to O(ε 0 ) is then the surprisingly simple where a cotangent and tangent cancelled in the Π I (K)Π J (P ) 2 term. Here, one can analyti- which only leaves the integrals in Ω to be computed numerically.
Let us now return to the χ ≈ 0 term in Eq. (97). We want to use tensor reduction to write theP integral inside of Ω in a simpler way. As written, theP integral in Eq. (97) depends on w =K ·P , so we cannot use full d-dimensional symmetry, but only a restricted (d−1)-dimensional symmetry, as detailed in Ref. [38]. However, we opt for another approach, which will yield a much simpler final expression.
Let us consider theP integral first in Eq. (97), so that we focus on Now expand out the denominator in a geometric series (here, we omit the Π(K) d−2 µν for space), resulting in Here, in going to from the first to second line, we used the symmetry of theP integral to remove all the terms with an odd number ofP . Thus, we must evaluate a generic integral of the form In Ref. [38], they provide the result 13 (translated to our notation) along with a recurrence relation between totally symmetric tensors 14 . Using the recurrence relation, and the Ward identity for Π(K), we conclude Additionally using where B is the Euler Beta function, we find in our case 13 see Eq. (30) in Ref. [38] 14 see Eq. (21) in Ref. [38] Therefore, Eq. (97) can be reduced to and Combining all of these results as in Eq. (90), we define Tr Π(P ) so that the total contribution from I A -I D , valid up to O(ε 0 ), reads: with ∆I A as in Eq. (103).

Extracting the coefficients of the ε expansion
We now turn to performing the ε expansions of [I ABCD (d)] UV above to calculate the terms in the expansion of the pressure in Eq. (32). The only divergence is contained in the overall trigonometric coefficient, which has the expansion π 2 sec πd 2 and so we see that Eq. (115) above contributes to p −2 in the expansion of the pressure. Moreover, from the ε expansion of the remainder of the integrand, this term likewise contributes to the remaining p −1 and p 0 terms as well. Recall that the Πs appearing in this expression are d-dimensional self energies, and thus to calculate the subleading contributions, we must expand not only the explicit d = 3 − 2ε appearing in the expression and the measure, but we also must expand the d-dimensional HTL self energies in a series The explicit details are shown in Appendix B 1. Notice that the ε −2 -divergence in Eq. (117) means that we will need the expansion up to Π I,2 in the calculation.
However, it turns out that we can avoid the explicit appearance of Π I,2 in our expressions, using the following approach. We write down a simpler integral [I −2 (d)] UV that will (i) remains. Using this as motivation, we define As shown in Appendix D, this integral can be performed analytically, yielding where ψ is the digamma function With these new definitions, we can thus split I ABCD into three pieces, which contribute to Eq. (32) in the following manner: Performing the ε expansion in the three terms identified in this list, we arrive at a set of angular integrals to perform, only some of which can be performed analytically. A full list of these contributing integrals is given in Appendix D. In Tab. II, we summarize the computed contributions to the coefficients p −2 , p −1 , and p 0 . As most of the one-dimensional integrals contributing to p 0 must be performed numerically, we mainly list numerical values for that row of the We now turn to the integral I F , given in Eq. (88), which we reproduce here: The evaluation of this integral proceeds in much the same way as the previous ones. However, we must first perform the contractions, which are more complicated in this case: as remarked above, the two 3g HTL vertices cannot be fully removed by the generalized Ward identies, since they are only contracted with resummed propagators.
Let us start by isolating the UV sensitivity in I F . It is straightforward to verify from the power countings in Tab. I and Eqs. (40) and (41) that I F has a UV logarithmic sensitivity only when all three momenta are hard. This leads to the following two conclusions for the UV sensitive part: Firstly, since the momenta K, P (and R) cannot be scale-separated when they are all hard, there is really only a single independent integration momentum (the variable X), and thus we expect only a single O(ε −1 ) divergence. Secondly, when all the momenta are hard, the resummed propagators approach their bare versions (see Eq. (39)), which leads to approximately direct contractions between the two 3g HTL vertices in this O(ε −1 ) term.
This second point in particular motivates us to isolate the "Kronecker-δ-like part" of the resummed propagators. We do this by unpacking the δ parts from the projection operators defined in Eqs. (B4)-(B5) and rearranging the terms to arrive at where I ∈ {T, L} and G TL (K) ≡ G T (K) − G L (K). Here, we have also introduced the with δ µµ T ≡ δ µi δ µ j δ ij and δ µµ L ≡ δ µ0 δ µ 0 . Hence, δ T picks out spatial indices and δ L picks out the temporal index. The UV-sensitivity will now arise from only taking the first term in Eq. (124) from each propagator. The two remaining terms are seen to be more supressed in the UV, and so lead to finite contributions.
Using the symmetry properties of the integrand KP R , we can now write I F in the form where the UV part [I F ] UV is given by The remaining finite parts ∆I (i) F are given by the following expressions

∆I
(1) and ∆I To obtain Eqs. (129) and (130), we first eliminated the 3g HTL vertices by using the generalized Ward identities and then plugged in the definition of G L from Eq. (27), where it occurs explicitly in Eq. (124). Note also that the term proportional to δΓ µνρ KP R K µ P ν R ρ has dropped out since it vanishes due to the Ward identity.
Before proceeding, we note the following fact, which we will use in evaluating these expressions. Namely, upon switching to the (X, χ) coordinates, the 3g HTL vertices scale where the δΓ µνρ KP R (χ, Ω) defined here depends only on the angles. This fact can be seen from the explicit integral representation of the 3g HTL vertex function in Appendix B 2.
We now proceed to the computation. As we did in the previous section, we shall only evalute the UV-sensitive term [I F ] UV in this section, and we defer the evaluation of the finite terms ∆I (i) F to Sec. III C below.

The X integral in d dimensions
We proceed to evalute [I F ] UV using the steps introduced at the beginning of this section.
After changing to the (X, χ) coordinates and using the scaling Eq. (131), we see that we have only one master integral in X to evaluate, namely where the function H d is given by the expression Here we have also introduced the notationŝ and the D-dimensional polar angle Φ R appearing in Π(R) is defined as Note that the divergence is contained in the trigonometric coefficient Using the above master integral, we find the following expression: Here, we have introduced a shorthand notation for contracted indices, e.g., [δΓ ijk KP R (χ, Ω)] 2 ≡ δΓ ijk KP R (χ, Ω)δΓ ijk KP R (χ, Ω). We note that in general d, the functions I d also depend on the angles χ and Ω.

Expanding in ε and doing the angular integrals
The details of the terms in this expansion are given in Appendix B 2, and the general approach is as follows: The 3g HTL vertex is given by a d-dimensional integral representation. We introduce a modified Feynman parametrization, which allows us to do these d-dimensional integrals, at least order-by-order in ε, leaving only the Feynman parameter left to be integrated over. Thus, we are able to write the functions in the above ε expansion as one-dimensional integral representations. As far as we are aware, such an explicit evaluation of the HTL vertices has never been performed before in the literature. Having discussed every truly UV-divergent integral, we are left with the finite terms.
These include the integrals I E and ∆I These terms all involve a coupling between the two loop-momenta to a degree that renders one unable to perform factorizations and other simplifications akin to those seen in Sec. III A. However, the finite nature of these terms makes them simpler to evaluate numerically, as there is no need to extract the divergent parts when evaluating them, unlike in the previously considered contributions. For this reason, we automate the evaluation of these terms, performing any remaining non-Lorentz-invariant tensor contractions by adapting the implementation discussed in Ref. [40] (although in Euclidean space).
Various properties of the HTL vertex functions, in particular the generalized Ward identities in Eq. (31) and their tracelessness, are again used extensively to reduce everything down to terms containing either no vertex corrections or the components δΓ 000 KP R or δΓ 0000 K,P,−K,−P , except in the case of ∆I (1) F . For this one separate term, we cannot fully remove the spatial components of the 3g HTL vertex function, but we reduce it as much as possible by usinĝ and then applying the generalized Ward identity on the first term. This way, we can reduce the number of spatial indices appearing in our expressions, as the 3g HTL vertices with more time components are easier to compute numerically. In the end, with this ∆I (1) F term, we can reduce all expressions to terms involving the four independent expressions δΓ 000 KP R ,  δΓ i00 KP R , δΓ ij0 KP R , and δΓ ijk KP R , with their indices contracted either with a second occurrence of one of these terms, or with the self energy functions Π 00 (Ŷ ), Π i0 (Ŷ ), and Π ij (Ŷ ), with Y ∈ {K, P, R}. More details about these manipulations can be found in Appendix C 1 c.
The 4g HTL vertex function appears in these finite terms for the first time; its properties are discussed in Appendix B 2. It should also be noted that, while here we are able to work in an integer dimension, computations of the HTL integrals in general d are easily automated using the method outlined in Appendix C 3, and in practice this method was applied in our computation.
With the integral scalarized, we simply follow the steps outlined in the beginning of Sec. III, setting d = 3 everywhere as these terms are finite. When changing to (X, χ) coordinates, we define new master integrals, which we do not show here but are easily computed. We also note here that similarly to the 3g HTL vertex, the 4g HTL vertex scales homogeneously as a function of X δΓ µνρσ K,P,−K,−P ≡ 1 X 2 δΓ µνρσ K,P,−K,−P (χ, Ω), which can be seen from the explicit integral form (Eq. (B36)) After this, the remaining nontrivial integrals (over χ, Φ K , Φ P , θ, and possible Feynman parameters in terms involving irreducible vertex corrections) are computed numerically using Monte Carlo methods [39].

IV. RESULTS AND CONCLUSIONS
We are now in a position to present our final results. By adding the corresponding elements of Tabs. II, III, and IV, we find the following fully soft contributions to the pressure in Eq. (32): p 0 = 17.150 (7).
Here, we note again for the reader that the angular integral Ω = π 2 /2 (see Eq. (80) for a full definition of the measure), and we have used a shorthand w = K · P/(|K||P |). The HTL self energy Π 0 and the 3g vertex correction δΓ 0 can be found in Appendices B 1 and C 1 respectively; note that we have also scaled out m E from them. In addition to these internal cross checks, we find that we reproduce the known coefficient of the leading logarithmic O(α 3 s ln 2 α s ) contribution to the pressure obtained in Ref. [20]. That work obtained this coefficient by expanding the two-loop HTL pressure in the semisoft region and using a cutoff prescription. Here, on the other hand, we used the fully resummed expressions, without expansions, and used dimensional regularization to arrive at the result.
Thus, the two techniques that yield agreement are quite independent. Note that if we use our result here for p −2 in the expansion in Eq. (33), we find that this result yields a coefficient for the O(α 3 s ln 2 α s ) term that is a factor of two larger than the result in Ref. [20] 15 . However, this is only the ln 2 (m E /Λ h ) term arising in p s 3 , which does not correspond to the full ln 2 α 1/2 s term in the sum over all regions. As we discuss in Appendix E using a simple example, it is natural to expect the ln 2 α 1/2 s coefficient to be precisely half of this ln 2 (m E /Λ h ) coefficient.
Thus, we do indeed find agreement between our result and Ref. [20].
As a further additional check with the literature, we have extended the formalism used in Ref. [20] to extract parts of the subleading logarithmic contribution, corresponding to pieces of p −1 . Within the cutoff regularization, there are two very distinct contributions to p −1 from the two-loop HTL pressure: One contribution arises from the soft-semisoft region, obtained by expanding the two-loop HTL pressure for only a single loop momentum, and corresponds to the second term of Eq. (142). The other contribution arises as a subleading correction to the semisoft-semisoft region (which in the (X, χ) coordinates has cutoffs on the radial X integral). This correction is regulator dependent, as it is sensitive to the exact ratio of cutoffs that occur in the double logarithm, rather than just the parametric size of the ratio. However, there is also a regulator-independent piece, corresponding in the language of the present paper to the last term of Eq. (142), due to the unique scaling properties of the doubly contracted vertex correction. As with the p −2 term, these regulator-independent pieces are found to be consistent between the two methods, giving us further confidence in our calculation.
We now turn to remarks upon our results. The first is that though we have rederived the same analytic value for the coefficient of the O(α 3 s ln 2 α s ) term, the integral expression for that result was much more complicated than the remarkably simple integral expression for p −2 above. In particular, in Ref. [20], the integral expression did not depend on only one momentum K, but rather both K and P . In fact, the integral given for p −2 above is exactly the same angular integral that appears in the Freedman-McLerran O(α 2 s ln α s ) result. Looking back through our analysis, we find that the angular tensor reduction was what lead to this vast simplification, as such a decomposition was not used in Ref. [20]. We further note that, since this same integral appears in p −2 , the corresponding contribution to p −1 [namely, the second term of Eq. (142)], is also an integral that appears in the α 2 s p s 2 term. 16 We also remark here upon the fact that in Ref. [20], the authors found that the same final result would be obtained if they set K on shell, with Π T,0 → Π T,0 (iK 0 = |k|, k) = 1/2 and Π L,0 → Π L,0 (iK 0 = |k|, k) = 0. This is also seen to work in the case of the O(α 2 s ln α s ) coefficient. (Here we have used our convention of factoring out m E .) However, we now see from the analysis in the present work that the replacement working for the O(α 3 s ln 2 α s ) coefficient follows directly from it working for the lower-order result, since the integral expressions are exactly the same. It is tempting to speculate that the leadinglogarithmic result at all orders may be related to the integral appearing in the α 2 s p s 2 term, and may allow for such a substitution.
As a further remark, we note that there have been many efforts to identify which diagrams give dominant contributions to an HTL calculation [42,43]. Here, we find, interestingly, that contributions containing irreducible HTL vertex corrections (i.e., those which cannot be removed using the generalized Ward identities) are clearly found to be smaller than those containing self energies. 17 We find that these corrections, though necessary for full, correct and by having the clear roadmap set forth in the Introduction, we believe that at least the further α 3 s p m 3 term can be obtained in a straightforward manner. In particular, there is no difficulty in combining the contributions from different kinematic regions, since they are all regulated consistently in dimensional regularization. With the α 3 s p m 3 term in hand, one could obtain the subleading O(α 3 s ln α s ) coefficient in the weak-coupling expansion, which may in turn allow one to use the principle of minimal sensitivity [44] to constrain the dependence of the pressure on the renormalization scaleΛ. This would potentially have important phenomenological implications for, e.g., the EOS of neutron-star matter [15,16]. Additionally, the entire organizational overview presented in the Introduction may have important consequences in itself: it may be possible to resum these logarithmic contributions to the pressure in some systematic way. An investigation of these points is left to later work.
Finally, we remark that there are some possible generalizations of this work which may be greatly aided by the organization and machinery that we have developed here. For example, including nonzero quark masses [45,46], or generalizations to nonzero temperature [19] might be possible using our present techniques. Such endeavors are, however, also left for the future. The standard group-theory factors for the SU(N c ) gauge group, some of which appear in the text, are given by where N c is the number of colors and f abc is the fully anti-symmetric structure constant. The generators of the fundamental representation (t a ) ij are normalized according to Tr t a t b = δ ab /2.
In the following we give the bare Feynman rules of QCD necessary for the computation of the diagrams under study. The vertices, with momentum flow towards the vertex, read: The corresponding propagators (with massless quarks), on the other hand, read We also define single-argument propagators to be two-argument propagators after enforcing the momentum-space delta function; for example (∆ 0 ) µν ab (P ) = δ µν δ ab P −2 − (1 − ξ)P µ P ν δ ab P −4 . In Sec. B below, we will write down some of the rules improved according to the standard HTL scheme.

Appendix B: The HTL framework
In this appendix, we outline in detail the Euclidean HTL framework used throughout this work. To make our HTL appendices self-contained, we repeat here some definitions that are scattered throughout the text. We leave the detailed evaluation of vertex functions to later appendices. We found the discussions in Refs. [25,37,47] helpful when creating this appendix.
is defined in the covariant gauge as where the parameter ξ fixes the gauge and The computation performed in this paper is most efficient to carry out in the ξ = 1 gauge, which we use throughout the text. The D-dimensional transverse and longitudinal projection operators, P µν T (K) and P µν L (K), are defined as with The Π T and Π L appearing in the denominators of the HTL propagators are components of the one-loop HTL self energy (or polarization tensor) Π µν ab ≡ δ ab Π µν of the gluon field. These functions satisfy the general relations and Π µµ (K) = (d − 1)Π T (K) + Π L (K), Additionally, the self energy satisfies the trivial Ward identity K µ Π µν (K) = 0.
b. Machinery for manipulating the propagator: We let ∆ µν 0 (K) denote the bare propagator with Π I (K) = 0, I ∈ {T, L}. We can extend this notation to label the other terms in the expansion of the full propagator in powers of the self energy: Here, we use a dot to represent contraction of adjacent indices, and we use the notation and make the identification [Π(K) 0 ] µν = δ µν to match the leading term. Thus, We can now introduce the following notation for the resummed propagator with the n leading terms removed: Consequently, the D µν n (K) are still resummed expressions, while the ∆ µν n (K) are not. Note that both ∆ µν n (K) and D µν n (K) are D-dimensionally transverse for every n, and that the following relations hold for any n ≥ 1: Notice also the full propagator at the end of both lines of Eq. (B14), and the fact that Eq. (B15) is not a partial sum of an infinite series, but is exact. The power-counting in Eq. (B17) allows one to use this notation to extract the form of the UV-sensitive terms in our calculations. c. Self energy: In the HTL approximation relevant for cold QM, the quark part of the one-loop gluon self energy is computed assuming that the momentum flowing along the quark lines is much larger than the external gluonic one. In this way, one obtains the result where we have introduced the lightlike four-vector V µ ≡ (−i,v) withv a unit vector in R d .

The integration measure in d dimensions is defined as
where z v ≡k ·v; note that the measure here is normalized to integrate to unity. The d-dimensional in-medium effective mass scale m E is given by This is the generalization of the effective mass scale to the case of multiple fermion flavors with different chemical potentials µ f at zero temperature. Throughout our text, m E denotes its d-dimensional value, and is never expanded in ε.
The scalar functions Π T and Π L can now be computed using the constraint equations in Eq. (B9) with the results where 2 F 1 is the hypergeometric function, and where the final equality assumes |k|/K 0 ∈ R and Re(d) > 1. If we now denote we find a very compact expression for the Π 00 integral with the notation Putting everything together, we find that the scalar functions Π T and Π L , expanded up to O(ε 2 ), can be expressed as where the coefficients above are given by and On occasion, we also denote Π n , n = 0, 1, 2, for the HTL self-energy tensor truncated to the appropriate order. It turns out to be convenient to express these results in terms of the polar angle Φ K , which is defined in Eq. (77). For example, in D = 4 dimensions, we obtain: Note that arctan[tan(Φ K )] = Φ K only for Φ K ∈ [0, π/2], Since we need expressions valid for the larger interval Φ K ∈ [0, π], we would need to replace arctan[tan(Φ K )] by Φ K − π · θ(Φ K − π/2), where θ denotes the Heaviside step function, if we wanted to further simplify Eq. (B28).

HTL effective vertices
As explained in the main text, treating the soft modes correctly within the HTL theory requires modifying not only the propagators, but also the n-point functions. This appendix contains the definitions of the three-(3g) and four-gluon (4g) vertices appearing at zero temperature.
a. The three-gluon vertex The effective 3g vertex is obtained by adding the HTL loop (which, at zero temperature, originates solely from the quark loop) to the bare vertex with the decomposition where the bare 3g vertex Γ µνρ 0 can be read off from Eq. (A2). The 3g HTL vertex function δΓ µνρ is in turn given by the expression The (tensor-valued) vertex function above is only defined when the sum of all of its arguments P, Q, and R is zero, and it is totally symmetric in its (Lorentz) indices (µ, ν, ρ) and traceless in any pair of indices, i.e. δ µν δΓ µνρ = 0 since V 2 = 0. Furthermore, it is even (odd) under even (odd) permutations of (P, Q,R).
Contracting Eq. (B31) with one of the momenta, for example with P µ , yields Comparing this to Eq. (B18), we find that the 3g HTL vertex function obeys the generalized Ward identity given in Eq. (31).

b. The four-gluon vertex
The effective 4g vertex is given by the decomposition where the bare 4g vertex (Γ 0 ) µνρσ abcd can be found in Eq. (A3). The general expression for the 4g HTL vertex function δΓ µνρσ abcd is uniquely determined from the knowledge of its symmetries, the 3g HTL vertex, as well as the Ward identities in Eq. (31). Here, however, we limit our detailed discussion only to the special case that we need in Eq. (22). This is, we take In this special case, the 4g HTL vertex δΓ µνρσ is given by the expression Akin to the 3g vertex correction, the 4g vertex function is totally symmetric in its four

Evaluation of the 3g HTL vertex function
We start the explicit evaluations of the HTL structures by considering the 3g vertex function δΓ µνρ P QR . The generalized Ward identities [see Eq. (31)] can often be used together with the tracelessness of the vertex correction to significantly simplify contributions containing the vertex correction. However, even for the two-loop HTL diagrams, the full structure of the HTLcorrected vertex is required due to the sunset diagram with two vertex corrections, as seen in Eq. (21). As a specific example, the contraction of two vertices (δΓ µνρ ) 2 ≡ δΓ µνρ δΓ µνρ18 includes every term allowed by the remnant d-dimensional rotational symmetry. As such, we must compute the following four independent vertex contributions δΓ 000 P QR , δΓ i00 P QR , δΓ ij0 P QR and δΓ ijk P QR . It turns out to be convenient to rewrite the expression in Eq. (B31) in the more symmetric form To evaluate the integral over the angles, one could try to combine the products in the two denominators into a single expression by using the "standard" Feynman parameterization However, in order to avoid the complications related to uQ 0 + (1 − u)(P 0 + Q 0 ) or −uQ 0 + (1 − u)(P 0 + Q 0 ) changing its sign at some value of u within the unit interval, causing the denominator of the integrand of Eq. (C2) to potentially vanish for some P and Q, we need to generalize the way that Feynman parameters are introduced. Let us first introduce the "symmetric" form of the parametrization to reach the general form where we have defined Evidently, the definition requires nonzero imaginary parts of Q·V and (P +Q)·V . 19 With this assumption, the modified parametrization can be shown to be equivalent to the standard form not only whenever the denominator of the latter is strictly nonvanishing, but also to yield 1/(AB) when the standard form does display divergent behaviour. We show this explicitly in Appendix C 4. Similarly, for the second term in Eq. (C1), we find The 3g HTL vertex function now takes the form where the four-vectors T and S are defined as: For further discussion, we will need the angular integral computed in Eq. (C44); we list here in d = 3 − 2ε dimensions some special cases as master integrals that will be used repeatedly in the following section where the master integrals are expanded to order ε and we have introduced the compact notation a(S) ≡ 1 + S 0 iL(S), b(S) ≡ ln(2) − 1 + ln (|s|/S) L(S) + 1 4|s| Li Here, the function L(S) is defined in Eq. (B24) and Li 2 is the dilogarithm function. 19 However, it can be generalized to momenta with vanishing zero-components; see Appendix C 4.
In the following subsections, we describe how to further evaluate the four independent vertex contributions δΓ 000 P QR , δΓ i00 P QR , δΓ ij0 P QR and δΓ ijk P QR . We also compute the O(ε) corrections to these vertex functions, which are needed in order to obtain the full O(ε 0 ) contribution to the UV term presented in Sec. III B.
a. The δΓ 000 P QR function Let us first concentrate on the δΓ 000 P QR vertex function. By using the general expression in Eq. (C6), we easily obtain where in the second line the v integral is performed by using a master integral listed in Eq. (C8). Correspondingly, the δΓ 000 P QR function squared can now be easily computed by using the expression above with two Feynman parameters u 1 and u 2 , where the variables T i and S i for i = 1, 2 are defined as follows: S i ≡ (1 + u i )P + (1 − u i )σ (P,P +Q) (P + Q).
(C12) b. The δΓ i00 P QR , δΓ ij0 P QR and δΓ ijk P QR functions We then proceed to describe the evaluation of the vertex functions δΓ i00 P QR , δΓ ij0 P QR and δΓ ijk P QR . By using the general expression in Eq. (C6), we obtain: The angular integrals appearing in Eq. (C13) can be dealt with by using the d-dimensional rotational symmetry. For instance, the rank-one integral can be written as Contracting both sides with the vector s i , and noting that we find for the reduction coefficient f 0 the following form After performing the remaining angular integrals in Eq. (C16) by using the master integrals listed in Eq. (C8), we find the result All in all, the vertex function δΓ i00 P QR then takes the form To evaluate the vertex functions δΓ ij0 P QR and δΓ ijk P QR further, we use the following tensorintegral reduction vv ivj (S · V ) 2 = δ ij f 00 (S 0 , |s|) + s i s j f 12 (S 0 , |s|), where the notation {δs} ijk ≡ δ ij s k + δ ik s j + δ jk s i has been introduced. Contracting Eq. (C19) with the Kronecker delta and vector s i , it is straightforward to show that the reduction coefficients above can be written as and f 000 (S) = 1 2|s| 4 3S 0 + (S 2 + 2S 2 0 )iL(S) Finally, inserting these results into Eq. (C13), we obtain for the vertex function δΓ ij0 P QR the following expression − P 0 σ (P,P +Q) δ ij f 00 (S) + s i s j f 12 (S) .
(C28) c. 3g HTL vertices contracted with external momenta In this section, we show how to evaluate the 3g HTL vertices contracted with external momenta. These techniques are used extensively in Sec. III C. Let us first consider the case where the 3g HTL vertex δΓ µνρ P QR is contracted with P µ T = δ µi p i . For this term, we cannot fully remove the spatial components of the 3g HTL vertex, but we can reduce it as much as possible using the identity To reduce the number of spatial indices appearing in our expressions above, we use the generalized Ward identity in Eq. (31) on the first term. This gives This method can be easily applied to the more complicated cases P µ T Q ν T δΓ µνρ P QR and P µ T Q ν T R ρ T δΓ µνρ P QR . By using Eq. (C29) and Ward identities, the following relations can be derived: and Here, the different components Π 00 (Y ), Π 0i (Y ) and Π ij (Y ) of the self energy Π µν (Y ) are given where Y ∈ {P, Q, R}. Note that the 3g HTL vertices with more time components are easier to compute numerically.

Evaluation of the 4g HTL vertex function
Following the discussion on the 3g HTL vertex correction, we will next consider the 4g vertex correction. A priori, it is considerably more complicated, and in order to handle the vertex in its full generality, a sensible option would be to turn to automation (see Appendix C 3). However, for the N 3 LO pressure, there is only a single resummed graph involving the 4g vertex correction, and it includes only a single vertex. It is easy to see that applying the symmetries and the Ward identities of the vertex correction along the same lines as in the previous section can only lead to a single irreducible term containing the 4g vertex correction δΓ 0000 P,Q,−P,−Q 20 : Following the 3g computation, we will apply a Feynman parametrization to make the numerics more tractable. However, as before, we must generalize the parametrization, following the discussion of Appendix C 4 [see also Eq. (C3)]. We combine the two common factors in the denominators via and to include the third factor, we denote to obtain 1 Hence, the full generalized Feynman parametrization for three-term denominator reads and analogously for P · V . Given this, δΓ 0000 P,Q,−P,−Q admits a representation δΓ 0000 P,Q,−P,−Q = 2 4 σ P +Q,P −Q where thev-integral can be obtained from Eq. (C44), and reads v 1 Substituting the leading-order term, we get δΓ 0000 P,Q,−P,−Q = 16σ P +Q,P −Q For us, setting ε = 0 suffices, as the 4g vertex correction only appears in the finite term Eq. (57).

Evaluating higher-rank HTL integrals
In Appendices C 1 and C 2 we have discussed special cases of calculations involving HTL vertices. However, for example for the purposes of automation and possible future, more complicated, computations, it is useful to be able to discuss the integrals that arise on a more general level. The prototypical tensor integrals arising in HTL calculations are of the  A, B). In the present paper, we have encountered the need to combine a factorized denominator to obtain a Feynman-like parametrization for integrals such as Eq. (C1). They involve arbitrary points A, B for which a nonvanishing denominator is not guaranteed.
Here we show that a suitable generalized Feynman parametrization is by means of a detailed proof that F (A, B) = 1/ (AB) for any A, B ∈ C\R. The missing case of real A, B will be briefly covered near the end of the section. In the following, we denote σ ≡ sgn(Im A Im B) in analogy with the σ defined in Appendix C 1.
First, we assume σ = +1. Then necessarily 0 / ∈ γ (A, B), so that we merely check that the identity holds in this standard case. Now, replace u → 1 − 2t, so that (−1, 1) → (0, 1) with the orientation reversed, leading to where we see a posteriori that the change of variables is permitted by assumption of 0 / ∈ γ (A, B). Note that this includes the special case B = A.
Next, consider σ = −1, this time without restrictions on whether or not origin is contained in γ (A, B), but by first assuming A + B = 0. This case requires slightly more care. For σ = −1 to be true, A, B must lie on opposite half-planes. As we exclude reals, barring the aforementioned A = −B we will have then covered all cases where 0 ∈ γ (A, B). As a consequence of the assumptions, the denominator must always have a nonzero imaginary part for all u ∈ (−1, 1), keeping it from vanishing on the interval. To see this, recall first that by assumption Im A = 0 (and Im B = 0). Should the denominator vanish for some u * ∈ (−1, 1), we would be lead to the equality The condition σ = −1 sets Im B/ Im A < 0, and we immediately see that there does not exist a u * ∈ (−1, 1) such that (1 + u * ) / (1 − u * ) ≤ 0, a contradiction. Therefore the integrand is where Θ(A) = Im A for Im A = 0 and Θ(A) = A otherwise. We will not cover this in detail, as it is unnecessary for us and the proof mostly consists of reapplying the above steps together with the implications σ = −1 leads to in these special cases. Furthermore, an extension to multiple denominators is straightforward by differentiation and an iterative application of the two-point formula, with an explicit example covered in Appendix C 2. In such cases, the geometric interpretation of the origin lying in the line segment connecting the two factors of the denominator is naturally generalised to the origin being within the convex hull of the points of the factorised denominator [48]. This condition serves as a check to see whether or not a generalized parametrization is necessary.
where Ω is as in Eq. (80). We define f (K, P ) 3 similarly in the obvious way. Observe that if f only depends on some proper subset of the variables integrated over, this still corresponds to the angular average of f over the angles upon which it depends.
We also note that in this appendix we use the convention where we have rescaled all momenta by m E and then removed the tildes as in the later parts of our main text. With the integrals that we will evaluate here, this is equivalent to simply setting m E = 1.

The integral in [I −2 (d)] UV
In this section, we demonstrate that which was used in Eq. (121) above. To perform this averaging, we use the definition of the self energy given in Eq. (B18), leading to where we use the lightlike four-vectors V defined in Appendix B 1 as well as the analogously defined U µ ≡ (−i,û) withû ∈ R d a unit vector. Additionally, we have defined the unit vectorN in the temporal direction. In going from the second to the third line above, we have changed variables in one term,û ↔v, to combine two terms. Now observe the identity where we used the definition of U , and the fact that i only appears in the temporal part.
But then inside an angular average overK, by multiplying the numerator and denominator byK · U * , we find We note that the function we are averaging over depends only on the components ofK within span(N ,û) (and we note that {N ,û} form a perpendicular basis for this subspace). If we splitK = K || + K ⊥ , with K || ∈ span(N ,û) and K ⊥ ·N = K ⊥ ·û = 0, then where we have recognized the denominator simply as |K || | 2 . We therefore see that the angular average involved above is the average of a unit vector over all of its directions, which simply leads to a constant multiple of the identity within its span δ µν || . However, U only depends on vectors within that span as well, and so we see Note that the entire analysis of this term could be conducted in the subspace span(N ,û).
Using the above results, our original average simplifies to We shall now analyze this within the three-dimensional subspace span(N ,û,v). Let us first splitK = K || + K ⊥ as above, but this time with K || ∈ span(N ,û,v) and K ⊥ ·N = K ⊥ ·û = K ⊥ ·v = 0. Becauseû andv are perpendicular toN , we can set up a three-dimensional coordinate system as depicted in Fig. 7 to perform the integral. Note that because of the geometry within this subspace, only the angle betweenû andv is in the d-dimensional spatial subspace. The calculation proceeds as where, we used the fact that û = 1. This is the desired result.