<figure><img src="image_1.png" /></figure>

RECEIVED: January 9, 2024

REVISED: February 20, 2024

ACCEPTED: February 23, 2024

PUBLISHED: March 19, 2024

# The force-force correlator at the hard thermal scale of hot QCD

Jacopo Ghiglieri<sup>D,a</sup> Philipp Schicho<sup>D,b</sup> Niels Schlusser<sup>D,c</sup> and Eamonn Weitz<sup>D,a</sup>

E-mail: jacopo.ghiglieri@subatech.in2p3.fr, schicho@itp.uni-frankfurt.de, niels.schlusser@unibas.ch, eamonn.weitz@subatech.in2p3.fr

**ABSTRACT:** High-energy particles traversing the Quark-Gluon plasma experience modified (massive) dispersion, although their vacuum mass is negligible compared to the kinetic energy. Due to poor convergence of the perturbative series in the regime of soft loop momenta, a more precise determination of this effective mass is needed. This paper continues our investigation on the factorisation between strongly-coupled infrared classical and perturbative ultraviolet behavior. The former has been studied non-perturbatively within EQCD by determining a non-local operator on the lattice. By computing the temperature-scale contribution to the same operator in 4D QCD at next-to-leading order (NLO), we remove the ultraviolet divergence of the EQCD calculation with an opposite infrared divergence from the hard thermal scale. The result is a consistent, regulator-independent determination of the classical contribution where the emergence of new divergences signals sensitivities to new regions of phase space. We address the numerical impact of the classical and NLO thermal corrections on the convergence of the factorised approach and on the partial applicability of our results to calculations of transport coefficients.

**KEYWORDS:** Effective Field Theories, Finite Temperature or Finite Density, Quark-Gluon Plasma, Thermal Field Theory

ARXIV EPRINT: 2312.11731

<footnote><sup>a</sup>SUBATECH, Université de Nantes, IMT Atlantique, IN2P3/CNRS, 4 rue Alfred Kastler, La Chantrerie BP 20722, 44307 Nantes, France</footnote>

<footnote><sup>b</sup>Institute for Theoretical Physics, Goethe Universität Frankfurt, 60438 Frankfurt, Germany</footnote>

<footnote><sup>c</sup>Biozentrum, University of Basel, Spitalstrasse 41, 4056 Basel, Switzerland</footnote>
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Background</td><td>4</td></tr><tr><td>2.1 Power counting and hierarchy of scales in jet-medium interactions</td><td>5</td></tr><tr><td>2.2 EQCD setup</td><td>5</td></tr><tr><td>3 Thermal modes at O(g<sup>2</sup>)</td><td>7</td></tr><tr><td>3.1 Equations of motion of the light-like Wilson line</td><td>8</td></tr><tr><td>3.2 Decomposition of the gluon self-energy</td><td>10</td></tr><tr><td>3.3 Z<sub>g,q<sup>-2</sup></sub><sup>(c)</sup> contribution</td><td>11</td></tr><tr><td>3.4 Z<sub>g,L-T</sub><sup>(c)</sup> contribution</td><td>12</td></tr><tr><td>4 Cancellation of IR QCD divergences with UV EQCD divergences</td><td>15</td></tr><tr><td>5 Summary and discussion</td><td>17</td></tr><tr><td>A Conventions</td><td>20</td></tr><tr><td>B One-loop self-energies and projectors</td><td>21</td></tr><tr><td>C Diagram (c) computation: technical details</td><td>22</td></tr><tr><td>C.1 Euclidean evaluation of Z<sub>g,q<sup>-2</sup></sub></td><td>22</td></tr><tr><td>C.2 Real-time evaluation of Z<sub>g,L-T</sub></td><td>27</td></tr><tr><td>D Determination of the QCD coupling</td><td>33</td></tr></table>

## 1 Introduction

The Quark-Gluon Plasma (QGP) is a deconfined, chirally symmetric phase of nuclear matter as described by Quantum Chromodynamics (QCD). The QGP is actively investigated through heavy-ion collisions. Disentangling the properties of this short-lived state from the shower of particles in the detectors requires an interplay of observables in the two main classes of bulk properties and hard probes. The former pertain to the collective properties of the many lower-energy particles, while the latter refer to the few energetic and/or feebly-coupled ones which are not expected to equilibrate.

Jets are one of the most important hard probes — see [1–3] for recent reviews. These self-collimated showers of hadrons are seeded by highly-energetic partons, in a regime where perturbative QCD should be applicable. Whether that remains true as these partons cross the QGP is still an open question, thus influencing the methods that one should apply to describe jets in medium. If the jet and the medium are weakly coupled — which does not necessarily imply a weakly-coupled medium — then the main theory ingredient in jet modification in
the QGP is medium-induced radiation. These modifications are sourced by the jet partons undergoing frequent exchanges of transverse momentum with the medium through elastic scatterings, encoded in the transverse scattering kernel $C(b_{\perp})$. The jet partons can also undergo forward scattering, whereby they are first scattered off their original momentum state and then back into it, thus shifting the dispersion relation by an amount named $m_{\infty}$, the asymptotic mass.

The transverse scattering kernel $C(b_{\perp})$ and $m_{\infty}$ are the most important medium-dependent ingredients in the description of medium-induced radiation, see e.g. [4]. In this paper, we continue our investigations of $m_{\infty}$, started in [5, 6]. We recall that these mass shifts are written in terms of the gauge condensate $Z_g$ and the fermion condensate $Z_f$ [7]

$$m_{\infty,q}^2 = g^2 C_F (Z_g + Z_f), \quad m_{\infty,g}^2 = g^2 C_A Z_g + 2g^2 T_F N_F Z_f, \quad (1.1)$$

where $m_{\infty,q}^2$ applies for quarks and $m_{\infty,g}^2$ applies for gluons. Here $C_F = (N_c^2 - 1)/(2N_c)$ is the quadratic Casimir of the quark representation, $C_A = N_c$ is the adjoint Casimir, $N_F$ is the number of light (Dirac) quark species, and $T_F = \frac{1}{2}$.

Eq. (1.1) can be understood as arising from integrating out the energy scale of the jet $E \gg T$ and truncating at first order in $T/E$. The matching coefficients are determined at first order in $g$ — we will return to this issue later. If the scale of the hard parton is $E \gtrsim T$ rather than $E \gg T$, higher orders in the $T/E$ expansion become relevant. This also makes clear that ours is a distinct problem from the one recently tackled in [8–11], namely the two-loop and power corrections to Hard Thermal Loops and their relation to the asymptotic mass for quasiparticles obeying $T \gg p \gg gT$.

The condensates $Z_g$ and $Z_f$ are non-local and have a gauge-invariant definition in terms of the correlators in the Hard Thermal Loop (HTL) action [7, 12]

$$Z_f = \frac{1}{2d_R} \langle \bar{\psi} \frac{v_\mu \gamma^\mu}{v \cdot D} \psi \rangle, \quad (1.2)$$

$$Z_g = \frac{1}{d_A} \langle v_\alpha F^{\alpha\mu} \frac{1}{(v \cdot D)^2} v_\nu F^\nu_\mu \rangle, \quad (1.3)$$

where $v^\mu = (1, \mathbf{v})$ is the light-like four-velocity of the hard particle, $d_R$ is the dimension of the fermion, $d_A$ the dimension of the adjoint representation, and the expectation value $\langle \dots \rangle$ denotes a thermal expectation value. Appendix A specifies our conventions.

All operators in thermal QCD are sensitive not just to the contribution of thermal modes, with momenta of the order of the temperature $T$. There exist also *soft modes*, with momenta $p \sim gT$, and gluonic ultrasoft modes, with momenta of $O(g^2T)$. The former cause the perturbative expansion to be in $g$ rather than $\alpha_s$ and are largely responsible for its slow convergence. The latter, at an operator-dependent order, cause the loop expansion to fail [13]. Moreover, for bosons these infrared (IR) modes are dominated by classical field dynamics, as their occupation numbers are $n_B(p) \approx T/p \gg 1$, with $n_B$ the Bose-Einstein distribution.

Soft modes first contribute to $Z_g$ at relative $O(g)$, whereas $Z_f$ is first expected to receive soft corrections at higher orders [7], i.e.

$$Z_g = Z_g^{\text{LO}} + \delta Z_g = \frac{T^2}{6} - \frac{T m_D}{2\pi} + \mathcal{O}(g^2), \quad Z_f = Z_f^{\text{LO}} + \delta Z_f = \frac{T^2}{12} + \mathcal{O}(g^2), \quad (1.4)$$
where $m_D^2$ is the leading-order Debye mass, given by [14]

$$m_D^2 = \frac{g^2 T^2}{3} (C_A + T_F N_F) + \mathcal{O}(g^4). \qquad (1.5)$$

The next-to-leading order (NLO) correction, $\delta Z_g$, is negative and, for values of the coupling corresponding to QGP temperatures, of comparable magnitude to the leading-order (LO) contribution $Z_g^{LO}$. The main motivation for [5, 6] and for the present paper is then to use an interplay of lattice and perturbation theory to determine the classical IR contribution to all orders. This can also allow to extend the applicability of the evaluation of the mass shift to lower temperatures, since it uses perturbation theory for the thermal modes only, which still have $\alpha_s$ as expansion parameters.

At the technical level, this has been made possible by developments introduced by Caron-Huot in [7, 15]. These papers showed how, for a hard, light-like parton interacting with the QGP, this classical contribution can be treated in Electrostatic QCD (EQCD) [16–20], a dimensionally-reduced Effective Field Theory (EFT) of QCD. This allows for non-perturbative evaluations of the EQCD contribution to $C(b_\perp)$ and $Z_g$ through lattice EQCD. We refer to [21–23] for the evaluation of $C(b_\perp)$ in EQCD, its merging with the T-scale contribution and its impact on medium-induced radiation.

For $m_\infty$, one can rely on the equivalence between inverse covariant derivatives and integrals over Wilson lines to rewrite eq. (1.3) as an integral over the light-like separation of the two field strengths, inserted over light-like Wilson lines — see eq. (2.1). In [5, 6] its EQCD counterpart was measured on the lattice as a function of the separation. To perform the separation integration, ref. [6] also provided a NLO perturbative determination in EQCD, to be used at distances shorter than the lattice spacing. There, in addition to discretisation, one faces another conceptual obstacle: the ultraviolet (UV) of EQCD is not, by construction, the UV of 4D thermal QCD. This is made manifest by the emergence of power-law and log-divergences from the integration of the LO (one-loop) and NLO (two-loop) perturbative EQCD at short separations. The super-renormalizable nature of EQCD prevents further UV divergences at higher orders.

In [6] the separation integral could then be carried out by subtracting the LO and NLO UV-divergent behaviors. For the power-law divergent LO term, this only amounts to avoiding double counting with the soft limit of the T-scale LO contribution $Z_g^{LO}$. For the log-divergent NLO term, this subtraction introduces a cutoff scale $L_{min}$. Physically, it signals an expected IR logarithmic divergence from thermal modes in 4D QCD at two-loop order.

The goal of this paper is precisely to determine the contribution to $Z_g$ from thermal modes in QCD at two-loop level. In so doing, we will recover the expected IR divergence, evaluate it in *dimensional regularisation* (DR) together with the aforementioned subtraction, thus lifting the dependence on the cutoff scale $L_{min}$ and providing an expression for $Z_g$ that is two-loop accurate for the T-scale modes and, when combined with the lattice results from [6], all-order accurate for the classical IR modes. We anticipate that we shall also find new ultraviolet and collinear divergences, signaling the emergence, at this loop order, of new regions of phase space and a potential sensitivity to one-loop corrections to the matching coefficients in eq. (1.1). Part of our main results have been anticipated in [24, 25]. We also address which part of our results here and in [5, 6] remains valid when the energy of the
<figure><img src="image_2.png" /><figcaption>Figure 1. Diagrams contributing the QCD force-force correlator $Z_g$ (2.1) at leading and next-to-leading order. An external gray shaded vertex denotes an $F^{-\perp}$ insertion; internal 2-point blobs the respective self-energy; a solid line a Wilson line; and a curly line a gauge boson ($A_\mu$).</figcaption></figure>

parton becomes of the order of the temperature, which is the relevant regime for calculations of transport coefficients and production rates.

The paper is organised as follows. Section 2 summarises the setup of the problem, the EQCD evaluation and the subtraction procedure for the UV divergences introduced in [6]. In section 3 we present our two-loop calculation of the T-scale contribution in 4D QCD, showing how, in addition to expected IR and vacuum UV divergences, we find extra ones and discuss their origin. In section 4, we show how the IR divergences cancel against the UV ones, leading to a finite, regulator-independent result for the classical, non-perturbative contribution. In section 5, we summarise our main findings and draw our conclusions. The appendices are dedicated to the conventions and technical details of our calculations.

## 2 Background

As we anticipated, $Z_g$ can be rewritten as the following integral over light-like separated field-strength insertions on a light-like Wilson line

$$Z_g = \frac{1}{d_R C_R} \int_0^\infty dx^+ x^+ \mathrm{Tr} \langle U_R(-\infty; x^+) v_\mu F^{\mu\nu}(x^+) U_R(x^+; 0) v_\rho F_\nu^\rho(0) U_R(0; \infty) \rangle, \quad (2.1)$$

where the Wilson line itself reads

$$U_R(x^+; 0) = \mathrm{P} \exp \left( ig \int_0^{x^+} dy^+ A^{-a}(y^+) T_R^a \right), \quad (2.2)$$

after having fixed $v^\mu = (1, 0, 0, 1)$ — see appendix A. The Wilson lines stretch from and back to negative light-cone infinity. They descend from the relation to the density matrix of the hard parton, which is defined in the asymptotic past.

At leading order, the Wilson lines source no gluons and one just connects the two field-strength tensors with a free thermal gluon propagator, as in diagram (a) of figure 1. It yields [26–28]

$$ \begin{align} Z_g^{\mathrm{LO}} &= (D-2) \int_Q [\theta(q^0) + n_B(|q^0|)] 2\pi\delta(Q^2) = (D-2) \int_Q n_B(|q^0|) 2\pi\delta(Q^2) \nonumber \\ &\stackrel{D=4}{=} 2 \int_q \frac{n_B(q)}{q} = \frac{T^2}{6}, \tag{2.3} \end{align} $$

where $n_B(q^0) \equiv (\exp(q^0/T) - 1)^{-1}$ is the Bose-Einstein distribution. Dimensional regularisation in $D = 4-2\epsilon$ ($d = 3-2\epsilon$ for the spatial dimensions) has been used to remove the scale-free vacuum contribution proportional to $\theta(q^0)$. See appendix A and eq. (A.2) for our conventions.
## 2.1 Power counting and hierarchy of scales in jet-medium interactions

The LO contribution to $Z_g$ arises from the scale $T$ as seen in eq. (2.3). In the infrared $n_B(q) \approx T/q$, so that eq. (2.3) has an order-$g$ sensitivity to the soft scale $gT$, viz. $\int_{q \sim gT} \frac{T}{q^2} \sim gT^2$. Indeed, that is where the NLO correction in eq. (1.4) appears. We can then summarise the perturbative expansion of $Z_g$ as follows [6]

$$Z_g = \left[ \frac{T^2}{6} - \frac{T\mu_h}{\pi^2} \right] \quad (2.4)$$

$$+ \left[ -\frac{Tm_D}{2\pi} + \frac{T\mu_h}{\pi^2} \right] \quad (2.5)$$

$$+ \left[ c_{\text{hard}}^{\ln} \ln \frac{T}{\mu_h} + c_T + c_{\text{hard}}^{\ln} \ln \frac{\mu_h}{m_D} + c_{\text{soft}}^{\ln} \ln \frac{m_D}{\mu_s} + c_{gT} + c_{\text{soft}}^{\ln} \ln \frac{\mu_s}{g^2 T} + c_{gT^2} \right] \quad (2.6)$$

$$+ \mathcal{O}(g^3). \quad (2.7)$$

We highlight the contributions from different orders of the coupling $g$ row by row and from different scales column by column. The $T$-scale contribution at $\mathcal{O}(g^0)$ from eq. (2.3) is displayed in (2.4) and the first $\mathcal{O}(g)$ correction from soft modes in (2.5). At the next order, all scales contribute: the completely non-perturbative color-magnetic scale $g^2T$, as well as the soft modes of the color-electric scale $gT$, and the thermal modes from $T$. The coefficients $c_{\text{hard}}^{\ln}$ and $c_{\text{soft}}^{\ln}$ were determined as a byproduct of [6], whose main aim was the non-perturbative, all-order determination of the classical soft and ultrasoft contribution. As we anticipated, the procedure used there introduced an artificial dependence on a cutoff $L_{\text{min}}$: for the moment we have kept the notation in eqs. (2.4)-(2.6) generic, with $\mu_h$ and $\mu_s$ taken as placeholders for an arbitrary scheme used to separate the thermal, soft and ultrasoft scales. The coefficients on eq. (2.6) are then to be understood as *scheme-dependent*. For reasons that will become evident later, knowing all coefficients in (2.4), (2.5), and (2.6) does not yield a self-contained estimate of $m_\infty$ at $\mathcal{O}(g^2)$, yet.

Consequently, our main aim is to connect the scheme used to remove the unphysical EQCD UV from the EQCD results in [6], giving rise to $L_{\text{min}}$ dependence, to dimensional regularisation. We will also compute the $T$-scale contribution in eq. (2.6), $c_{\text{hard}}^{\ln}$ and $c_T$, in that same scheme. To this end, we first briefly review the determination of the EQCD contribution to $Z_g$ from [5, 6].

## 2.2 EQCD setup

One starts from the EQCD continuum action,

$$S_{\text{EQCD}} = \int_{\mathbf{x}} \left\{ \frac{1}{2} \text{Tr} F_{ij} F_{ij} + \text{Tr} [D_i, \Phi][D_j, \Phi] + m_D^2 \text{Tr} \Phi^2 + \lambda_E (\text{Tr} \Phi^2)^2 \right\}, \quad (2.8)$$

where $\Phi$ is the adjoint, massive scalar arising from dimensional reduction of the $A^0$ gauge field. One then derives the non-unitary EQCD analog of eq. (2.1)

$$\tilde{U}_R(L; 0) = P \exp \left( ig_{3d} \int_0^L dz (A_z^a(z) + i\Phi^a(z)) T_R^a \right). \quad (2.9)$$
Since the field operators are classical in EQCD, they commute and allow us to write

$$ \frac{1}{d_R C_R} \mathrm{Tr} \langle \tilde{U}_R(-\infty; L) \mathcal{O}^a(L) T_R^a \tilde{U}_R(L; 0) \mathcal{O}^b(0) T_R^b \tilde{U}_R(0; -\infty) \rangle = \frac{T_R}{d_R C_R} \langle \mathcal{O}^a(L) \tilde{U}_A^{ab}(L; 0) \mathcal{O}^b(0) \rangle. \qquad (2.10) $$

Thus, the EQCD equivalent of $Z_g$ reads, given rotational invariance in the transverse plane

$$ Z_g^{3d} = -\frac{(d-1)T}{d_A} \int_0^\infty dL \, L \langle (F_{xz}^a(L) + i(D_x\Phi(L))^a) \tilde{U}_A^{ab}(L; 0) (F_{xz}^b(0) + i(D_x\Phi(0))^b) \rangle, \qquad (2.11) $$

where we used $d_R C_R / T_R = d_A$. At LO in perturbative EQCD, we can connect the operators with EQCD propagators to find [7]

$$ Z_{g LO}^{3d} = T \int_{\mathbf{q}} \frac{1}{(q^z + i\epsilon)^2} \left[ -\frac{q_{\perp}^2}{q^2 + m_D^2} + \frac{q_{\perp}^2}{q^2} + \frac{(d-1)q_z^2}{q^2} \right] = -\frac{Tm_D}{2\pi}, \qquad (2.12) $$

where the massive propagator is the $\Phi$ propagator. Here we used DR to get to a finite result: in so doing, a linearly-UV divergent term gets discarded. It is instructive to note that this term would be $(d-1) \int_{\mathbf{q}} \frac{T}{q^2}$, which corresponds precisely to the IR sensitivity of eq. (2.3), as argued in the previous subsection, and to the linear terms in $\mu_h$ in eqs. (2.4) and (2.5). Thus, dimensional regularisation corresponds to subtracting off that $(d-1) \int_{\mathbf{q}} \frac{T}{q^2}$ bare, soft limit from the resummed result in eq. (2.12). As we shall now see, this does not happen automatically in the scheme that follows from the lattice computation.

For that end, it is more convenient to set $d=3$ and to work with the quantities

$$ \langle EE \rangle = \frac{1}{2} \langle (D_x \Phi(L))^a \tilde{U}_A^{ab}(L, 0) (D_x \Phi(0))^b \rangle, \qquad (2.13) $$

$$ \langle BB \rangle = \frac{1}{2} \langle F_{xz}^a(L) \tilde{U}_A^{ab}(L; 0) F_{xz}^b(0) \rangle, \qquad (2.14) $$

$$ i \langle EB \rangle = \frac{i}{2} \langle (D_x \Phi(L))^a \tilde{U}_A^{ab}(L; 0) F_{xz}^b(0) \rangle + \frac{i}{2} \langle F_{xz}^a(L) \tilde{U}_A^{ab}(L; 0) (D_x \Phi(0))^b \rangle, \qquad (2.15) $$

which leaves us with

$$ Z_g^{3d} = -\frac{4T}{d_A} \int_0^\infty dL \, L \left( -\langle EE \rangle + \langle BB \rangle + i\langle EB \rangle \right) = \frac{4T}{d_A} \int_0^\infty dL \, L \langle FF \rangle, \qquad (2.16) $$

where we have implicitly defined $\langle FF \rangle = \langle EE \rangle - \langle BB \rangle - i\langle EB \rangle$.

In [5, 6] continuum-extrapolated lattice EQCD determinations of the three components of $\langle FF \rangle$ have been performed for distances $0.25 < g_{3d}^2 L < 3.0$ at four values of the temperature. To integrate over all separations, this lattice determination must be complemented at both long and short distances to provide a finite result for the all-order classical contributions from the soft and ultrasoft scales. In more detail, the IR region is addressed with a fitting ansatz, intended to reproduce the expected exponential falloff from electrostatic and magnetostatic screening, since one cannot reach arbitrarily large distances on the lattice. As for the UV, the lattice approach becomes impractical at short distances due to associated discretisation effects. However, at the shortest available separations, where $m_D L \ll 1$, the lattice determination agrees well with the perturbative NLO EQCD determination, as expected. This then suggests a switch to the perturbative EQCD determination at short distances. Yet, the
naive integration down to zero separation of the perturbative EQCD result would result in power-law and log divergences, as we mentioned.

But how does this divergent behavior manifest itself? Because of the super-renormalizability of EQCD, in the UV, it is natural to expect from dimensional analysis that

$$ \langle FF \rangle_{m_D L \ll 1} = \frac{c_0}{L^3} + \frac{c_2 g_{3d}^2}{L^2} + \frac{c_4 g_{3d}^4}{L} + \dots \qquad (2.17) $$

The first and second terms above are divergent when inserted into eq. (2.16). If integrated down to a $L_{min} > 0$, the former gives rise to a $c_0/L_{min}$ contribution. In this position-space scheme, this corresponds to the linear term in $\mu_h$ in eq. (2.5). As the previous discussion should have made clear, such a term corresponds to the soft limit of the bare LO calculation in eq. (2.3); as it is already included there, it should be subtracted off from the EQCD calculation, as per [5].

The coefficient $c_2$ was computed in [6] and, being a logarithmic divergence, shows up also in DR. To arrive at a determination for $Z_g^{3d}$, [6] merged together the perturbative, intermediate and IR regions as follows

$$ Z_g^{3d}|_{\text{merge}} = \frac{4T}{d_A} \left\{ \int_0^{L_{\text{min}}} dL \, L \left[ \langle FF \rangle_{\text{NLO}} - \frac{c_0}{L^3} - \frac{c_2 g_{3d}^2}{L^2} \right] + \int_{L_{\text{min}}}^{L_{\text{max}}} dL \, L \left[ \langle FF \rangle_{\text{lat}} - \frac{c_0}{L^3} \right] \right. \\ \left. + \int_{L_{\text{max}}}^{\infty} dL \, L \left[ \langle FF \rangle_{\text{tail}} - \frac{c_0}{L^3} \right] \right\}, \qquad (2.18) $$

where [5, 6]

$$ c_0 = \frac{d_A}{4\pi}, \qquad c_2 = \frac{C_A d_A}{(4\pi)^2}. \qquad (2.19) $$

The shortest and longest separations for which lattice data is available are $L_{min} = 0.25/g_{3d}^2$ and $L_{max} = 3.0/g_{3d}^2$.<fnref n="1" /> Between these two separations, the lattice-determined $\langle FF \rangle_{\text{lat}}$ is used. At asymptotically large distances the modeled IR tail $\langle FF \rangle_{\text{tail}}$ is used and at short distances the NLO perturbative EQCD result $\langle FF \rangle_{\text{NLO}}$ is integrated. As we mentioned, the $c_0$ term must be subtracted off everywhere, whereas the $c_2$-term is only subtracted in the perturbative region. Eq. (2.18) thus introduces an artificial logarithmic dependence on its boundary $L_{min}$.

In section 4 we will show how this subtraction can be recast in dimensional regularisation. The resulting UV pole cancels against an opposite IR one from computing the two-loop contribution to $Z_g$ in thermal 4D QCD, which shall be carried out in the next section. This is then our main motivation, as their sum provides a finite result for the non-perturbative contribution, with the artificial cutoff, $L_{min}$, lifted. Such a calculation is evidently also in line with our ultimate goal of determining all $O(g^2)$ corrections to $Z_g$.

### 3 Thermal modes at $O(g^2)$

We embark on the calculation of the T-scale contribution to $Z_g$ at NLO. In principle, this requires the computation of the two-loop diagrams (b)-(j) in figure 1. As we are dealing with

<footnote n="1">¹The $\langle EB \rangle$ contribution was only available for $L > 0.5/g_{3d}^2$ in [6]. This is irrelevant from the standpoint of $c_2$ subtraction, since only $\langle EE \rangle$ is logarithmically divergent at small $L$.</footnote>
the temperature scale, there is no need to or benefit in using any IR-resummed propagators and vertices in the real-time formalism. The unresummed nature of the propagators can be used to our advantage to greatly simplify the calculation.

In a nutshell, all but diagrams (c), (f), (g) only contain terms that are proportional to either $\langle A^- A^- \rangle$ or $\langle A^- A^x \rangle$, both of which vanish in Feynman gauge. Furthermore, following an argument to be given in section 3.1, we will see that only diagram (c) needs to be computed. The remainder of the section is then devoted to the splitting up and the computation of the different contributions: quark and gauge (gluon and ghost) in the polarisation tensor.

## 3.1 Equations of motion of the light-like Wilson line

Looking at eq. (2.1), we can write

$$Z_g = \frac{-1}{d_R C_R} \int_0^\infty dx^+ x^+ \text{Tr} \langle U_R(-\infty; x^+) F^{-\perp}(x^+) U_R(x^+; 0) F^{-\perp}(0) U_R(0; -\infty) \rangle, \quad (3.1)$$

where $F^{-\perp}F^{-\perp} = F^{-x}F^{-x} + F^{-y}F^{-y}$. Continuing onwards, let us recast eq. (3.1) as a double integral, using translation invariance

$$\begin{aligned} Z_g &= \frac{-1}{d_R C_R} \int_0^\infty dx^+ \int_0^{x^+} dy^+ \\ & \quad \times \text{Tr} \langle U_R(-\infty; x^+ + y^+) F^{-\perp}(x^+ + y^+) U_R(x^+ + y^+; y^+) F^{-\perp}(y^+) U_R(y^+; -\infty) \rangle, \\ &= \frac{-1}{d_R C_R} \int_0^\infty dx^+ \int_0^{\frac{x^+}{2}} dy^+ \text{Tr} \langle U_R(-\infty; x^+) F^{-\perp}(x^+) U_R(x^+; y^+) F^{-\perp}(y^+) U_R(y^+; -\infty) \rangle, \end{aligned} \quad (3.2)$$

where we have shifted $x^+$ in going to the second line. We now rewrite $Z_g$ as

$$Z_g \equiv Z_g^{\perp\perp} + Z_g^{\perp-} + Z_g^{-\perp} + Z_g^{--}, \quad (3.3)$$

with the index labelling to be elaborated on momentarily. By calling on a method introduced in [29], we rewrite

$$F^{-\perp} = -\partial^{\perp} A^{-} + [D^{-}, A^{\perp}]. \quad (3.4)$$

The equation of motion of the Wilson line, $D_{x^+}^{-} U(x^+; y^+) = 0$ can then be used to rewrite the commutator as a total derivative, i.e.

$$U_R(a; x^+) [D^-, A^\perp(x^+)] U_R(x^+; b) = \frac{d}{dx^+} [U_R(a; x^+) A^\perp(x^+) U_R(x^+; b)]. \quad (3.5)$$

This should shed light on the labelling in eq. (3.3): the indices there denote whether the $-\partial^\perp A^\perp$ ($\perp$) or the $[D^-, A^\perp]$ component ($-$) of the field strength tensor has been taken. We start with the first term from eq. (3.3),

$$Z_g^{\perp\perp} = - \frac{1}{d_R C_R} \int_0^\infty dx^+ x^+ \text{Tr} \langle U_R(-\infty; x^+) \partial^\perp A^\perp(x^+) U_R(x^+; 0) \partial^\perp A^\perp(0) U_R(0; -\infty) \rangle. \quad (3.6)$$
For the second term, we have, employing eq. (3.5)

$$Z_{g}^{\perp -} = \frac{1}{d_{R}C_{R}} \int_{0}^{\infty} dx^{+} \mathrm{Tr} \left\langle U_{R}(-\infty; x^{+}) \partial^{\perp} A^{-}(x^{+}) U_{R}(x^{+}; x^{+}/2) \right. \\ \left. \times \left[ A^{\perp}(x^{+}/2) U_{R}(x^{+}/2; 0) - U_{R}(x^{+}/2; 0) A^{\perp}(0) \right] U_{R}(0; -\infty) \right\rangle. \quad (3.7)$$

Translation invariance dictates that for the first part,

$$\begin{aligned} & \left\langle U_R(-\infty; x^+) \partial^\perp A^-(x^+) U_R(x^+; x^+/2) A^\perp(x^+/2) U_R(x^+/2; -\infty) \right\rangle \\ &= \left\langle U_R(-\infty; x^+/2) \partial^\perp A^-(x^+/2) U_R(x^+/2; 0) A^\perp(0) U_R(0; -\infty) \right\rangle, \end{aligned} \quad (3.8)$$

so that, by aptly redefining the integration variable, one finds, including its counterpart, that

$$Z_{g}^{\perp -} = \frac{1}{d_{R}C_{R}} \int_{0}^{\infty} dx^{+} \mathrm{Tr} \left\langle U_{R}(-\infty; x^{+}) \partial^{\perp} A^{-}(x^{+}) U_{R}(x^{+}; 0) A^{\perp}(0) U_{R}(0; -\infty) \right\rangle, \quad (3.9)$$

$$Z_{g}^{-\perp} = \frac{-1}{d_{R}C_{R}} \int_{0}^{\infty} dx^{+} \mathrm{Tr} \left\langle U_{R}(-\infty; x^{+}) A^{\perp}(x^{+}) U_{R}(x^{+}; 0) \partial^{\perp} A^{-}(0) U_{R}(0; -\infty) \right\rangle, \quad (3.10)$$

where we used the fact that, in a non-singular gauge, $A^{\perp}(\infty) = 0$ and we applied translation invariance and relabeled $y^{+}$ into $x^{+}$. Finally, for the term with both commutators, we find

$$Z_{g}^{--} = \frac{1}{d_{R}C_{R}} \mathrm{Tr} \left\langle U_{R}(-\infty; 0) A^{\perp}(0) A^{\perp}(0) U_{R}(0; -\infty) \right\rangle, \quad (3.11)$$

where we have again set the boundary term at infinity to zero.

At two-loop level, eq. (3.6) can only receive a contribution from diagram (c) in Feynman gauge; terms with an $A^{-}$ pulled from the Wilson line cannot contribute, as there is no triple $A^{-}$ vertex. Moreover, there cannot be a contribution of the form of diagram (g) as the $gA^{-}A^{\perp}$ part of the covariant derivative has been removed from the correlator through the previous manipulations. Similar arguments can be safely applied to eqs. (3.9) and (3.10): they only receive a “mixed” contribution from diagram (c); this relies on the observation that no vertex with two $A^{-}$ gluons exists.

Eq. (3.11) is the only piece of the decomposition that does not vanish at LO in Feynman gauge; it corresponds precisely to diagram (a). Eq. (3.11) also receives an NLO contribution from (c). It also contains the only contribution directly involving the Wilson lines, as either of them can source an $A^{-}$ gluon that can connect with a three-gluon vertex to the two $A^{\perp}$ ones. However, that would give rise to the following colour structure

$$C = f^{abc} A^{\perp a} A^{\perp b} A^{-c}, \quad (3.12)$$

which vanishes, due to the symmetric (anti-symmetric) nature of $A^{\perp a}A^{\perp b}$ ($f^{abc}$). Therefore, through this reorganisation, we conclude that all of the NLO corrections to $Z_{g}$ come from diagram (c). We arrive at the same conclusion from an explicit real-time computation of diagrams (f) and (g), whose contributions end up cancelling each other [25].
## 3.2 Decomposition of the gluon self-energy

We thus dive straight into the computation of diagram (c)

$$Z_g^{(c)} = -\frac{(D-2)}{d_R C_R} \int_0^\infty dx^+ x^+ \mathrm{Tr} \langle (\partial^- A^x(x^+) - \partial^x A^-(x^+))(\partial^- A^x(0) - \partial^x A^-(0)) \rangle, \quad (3.13)$$

where we have used rotational invariance in the transverse plane. As the field-strength insertions are at positive time separation, there is indistinguishably a forward Wightman or time-ordered one-loop propagator stretching between them

$$\begin{aligned} Z_g^{(c)} &= -(D-2) \int_0^\infty dx^+ \int_Q x^+ e^{-iq^- x^+} [(q^-)^2 G_>^{xx}(Q) - 2q^x q^- G_>^{x-}(Q) + (q^x)^2 G_>^{--}(Q)] \\ &= (D-2) \int_Q \frac{1+n_B(q^0)}{(q^- - i\epsilon)^2} (F_R(Q) - F_A(Q)), \end{aligned} \quad (3.14)$$

where we have used the Kubo-Martin-Schwinger (KMS) relation to write the forward Wightman > propagator in terms of the statistical function $1 + n_B(q^0)$ and the one-loop spectral function $F_R(Q) - F_A(Q)$. Its retarded and advanced components read

$$F_{R/A}(Q) = \frac{i}{(Q^2 \pm i\epsilon q^0)^2} \left[ (q^-)^2 \Pi_{R/A}^{xx}(Q) - 2q^x q^- \left( \Pi_{R/A}^{0x}(Q) - \Pi_{R/A}^{zx}(Q) \right) \right. \\ \left. \qquad + (q^x)^2 \left( \Pi_{R/A}^{00}(Q) + \Pi_{R/A}^{zz}(Q) - 2\Pi_{R/A}^{z0}(Q) \right) \right], \quad (3.15)$$

when written in terms of the one-loop self-energies, given in appendix B. It is then useful to decompose the self-energies in terms of their longitudinal (L) and transverse (T) components (with respect to the three-momentum), which can also be found in appendix B

$$\Pi_{\mu\nu}^{R/A}(Q) = P_{\mu\nu}^{L}(Q)\Pi_{L}^{R/A}(Q) + P_{\mu\nu}^{T}(Q)\Pi_{T}^{R/A}(Q) \quad (3.16)$$

where $P_{\mu\nu}^L$ and $P_{\mu\nu}^T$ are given in eqs. (B.1)–(B.2).

We can then rewrite eq. (3.15) as

$$F_{R/A}(Q) = \frac{i}{(Q^2 \pm i\epsilon q^0)^2} \left[ (q^-)^2 \Pi_T^{R/A}(Q) + \frac{(q^x)^2 Q^2}{q^2} \left( \Pi_L^{R/A}(Q) - \Pi_T^{R/A}(Q) \right) \right], \quad (3.17)$$

which leaves us with the decomposition

$$\begin{aligned} Z_g^{(c)} &= \int_Q \frac{1+n_B(q^0)}{(q^- - i\epsilon)^2} \left\{ \frac{i}{(Q^2+i\epsilon q^0)^2} \left[ (D-2)(q^-)^2 \Pi_T^R(Q) + \frac{q_\perp^2 Q^2}{q^2} \left( \Pi_L^R(Q) - \Pi_T^R(Q) \right) \right] - \mathrm{adv.} \right\} \\ &= Z_{g,q^{-2}}^{(c)} + Z_{g,L-T}^{(c)}, \end{aligned} \quad (3.18)$$

where adv. stands for the advanced part while $Z_{g,q^{-2}}^{(c)}$ corresponds to the $(q^-)^2$- and $Z_{g,L-T}^{(c)}$ to the $\Pi_L - \Pi_T$-proportional terms

$$Z_{g,q^{-2}}^{(c)} = \int_Q [1 + n_B(q^0)] \left\{ \frac{i(D-2)\Pi_T^R(Q)}{(Q^2 + i\epsilon q^0)^2} - \mathrm{adv.} \right\}, \quad (3.19)$$

$$Z_{g,L-T}^{(c)} = \int_Q \frac{1 + n_B(q^0)}{(q^- - i\epsilon)^2} \left\{ \frac{i q_\perp^2}{(Q^2 + i\epsilon q^0) q^2} \left( \Pi_L^R(Q) - \Pi_T^R(Q) \right) - \mathrm{adv.} \right\}. \quad (3.20)$$

The strategy adopted for the remainder of this section is to compute $O(g^2)$ corrections to $Z_g$ according to the above decomposition. Rather than wading through technical details (which can be found in appendix C), we instead focus on highlighting and discussing the physical relevance of the various divergences that we encounter.
### 3.3 $Z_{g,q-2}^{(c)}$ contribution

For the computation of $Z_{g,q-2}^{(c)}$, we can resort to standard Euclidean techniques. We start from

$$Z_{g,q-2}^{(c)} = i(D-2) \int_Q \left(\frac{1}{2} + n_B(q^0)\right) \left[ \frac{\Pi_T^R(Q)}{(Q^2 + i\epsilon q^0)^2} - \frac{\Pi_T^A(Q)}{(Q^2 - i\epsilon q^0)^2} \right], \quad (3.21)$$

where we have kept the even-in-$q^0$ part of the integrand only. By analyticity, we have

$$Z_{g,q-2}^{(c)} = -(D-2) \oint_{Q_E} \frac{\Pi_T^E(Q_E)}{Q_E^4}. \quad (3.22)$$

Here, $Q_E^2 = q_0^2 + q^2$ is the Euclidean four-momentum squared, with $q_0 = 2\pi nT$ a bosonic Matsubara frequency and $\oint_{Q_E}$ the standard sum-integral, see eq. (A.2). Likewise, $\Pi_T^E(Q_E)$ is the transverse Euclidean polarisation tensor, given in appendix B. The integrations can then be carried out using standard Euclidean techniques, as detailed in appendix C.1.

The contribution of a gluon and ghost loop reads (cf. eq. (C.18))

$$Z_{g,q^{-2},B,n \neq 0}^{(c)} = \frac{g^2 C_A T^2}{96 \pi^2} \left[ \frac{11}{3} \left( \frac{1}{\epsilon} + 4 \ln \frac{\bar{\mu}}{T} \right) + 8 \ln A - 4 \gamma_E + \frac{61}{9} + \frac{10}{3} \ln \pi - \frac{52}{3} \ln 2 + \frac{4}{3} \times [0.299378] \right], \quad (3.23)$$

where B stands for bosonic, $\gamma_E$ is the Euler-Mascheroni constant and $A \approx 1.28243$ is the Glaisher constant. Here and in all following formulae we truncate the expansion in $D = 4-2\epsilon$ to the order-$\epsilon^0$ term included. We have further subtracted off the $q_0 = 0$ zero-mode contribution associated with the $Q_E$ sum-integral. This IR-divergent zero-mode contribution can be found in section 4 (see eq. (4.1)), where we cancel its divergence with an opposite one from the EQCD subtraction term discussed in and after eq. (2.18).

The fermionic contribution from the quark loop reads — see eq. (C.19)

$$Z_{g,q^{-2},F}^{(c)} = -\frac{g^2 T_F N_F T^2}{96 \pi^2} \left[ \frac{4}{3} \left( \frac{1}{\epsilon} + 4 \ln \frac{\bar{\mu} A^3}{4\pi T} \right) + 4\gamma_E + \frac{26}{9} - \frac{10}{3} \times [0.199478] \right]. \quad (3.24)$$

In this case, the $Q_E$ zero mode is included, as its contribution is finite as expected from the purely non-Abelian nature of the $c_2$ subtraction in eqs. (2.18)–(2.19).

Upon summing eqs. (3.23) and (3.24), we find

$$\begin{aligned} Z_{g,q^{-2}}^{(c),\text{IR safe}} &= Z_{g,q^{-2},B,n \neq 0}^{(c)} + Z_{g,q^{-2},F}^{(c)} \\ &= \frac{g^2 T^2}{96 \pi^2} \left[ \left( \frac{11}{3} C_A - \frac{4}{3} T_F N_F \right) \left( \frac{1}{\epsilon} + 4 \ln \frac{\bar{\mu}}{T} \right) \right. \\ &\qquad \left. + C_A \left( 8 \ln A - 4 \gamma_E + \frac{61}{9} + \frac{10}{3} \ln \pi - \frac{52}{3} \ln 2 + \frac{4}{3} \times [0.299378] \right) \right. \\ &\qquad \left. - T_F N_F \left( 16 \ln A + 4 \gamma_E + \frac{26}{9} - \frac{16}{3} \ln(4\pi) - \frac{10}{3} \times [0.199478] \right) \right]. \end{aligned} \quad (3.25)$$

The attentive reader will have noticed that eqs. (3.23) and (3.24) are divergent. The very attentive reader will also have noticed that the divergent term in eq. (3.25) is proportional to
the first coefficient of the QCD $\beta$ function. Indeed, the divergence is of ultraviolet origin: it arises from the thermal part of a bosonic loop integral multiplying the vacuum part of the other loop integral. This divergence is then removed by charge renormalisation. In more detail, let us consider the standard relation between bare and renormalised coupling at one-loop level, viz.

$$g_{\text{bare}}^2 = g^2(\bar{\mu}) \left[ 1 - \frac{g^2 b_0}{(4\pi)^2 \epsilon} \right] + \mathcal{O}(g^6), \quad \text{where} \quad b_0 = \frac{11}{3} C_A - \frac{4}{3} T_F N_F. \qquad (3.26)$$

Let us then go back to the leading-order term in eq. (2.3), recalling that, as per eq. (1.1), $m_\infty^2$ is obtained by multiplying $Z_g$ by $g^2$.<fnref n="2" /> Hence, if that $g^2$ is replaced by eq. (3.26), the $\mathcal{O}(g^4)$ term there generates the counterterm

$$Z_{\text{g}}^{\text{LO,ct}} = -(d-1) \frac{g^2 b_0}{(4\pi)^2 \epsilon} \int_{\mathbf{q}} \frac{n_{\text{B}}(q)}{q} = \frac{g^2 T^2}{96\pi^2} \left( \frac{11}{3} C_A - \frac{4}{3} T_F N_F \right) \left[ -\frac{1}{\epsilon} - 2 \ln \frac{\bar{\mu} A^{12}}{4\pi T} + 1 \right], \quad (3.27)$$

when multiplying the LO term in eq. (2.3). Upon adding eqs. (3.25) and (3.27), we find<fnref n="3" />

$$\begin{aligned} Z_{g,q^{-2}}^{(c),\text{IR/UV safe}} &= Z_{g,q^{-2}}^{(c),\text{IR safe}} + Z_{g}^{\text{LO,ct}} \\ &= \frac{g^2 T^2}{96\pi^2} \left\{ \left( \frac{11}{3} C_A - \frac{4}{3} T_F N_F \right) 2 \ln \frac{\bar{\mu}_{\text{UV}}}{T} \right. \\ & \qquad \left. + C_A \left[ -80 \ln A - 4 \gamma_E + \frac{94}{9} + \frac{32}{3} \ln \pi - \frac{8}{3} \ln 2 + \frac{4}{3} \times [0.299378] \right] \right. \\ & \qquad \left. + T_F N_F \left[ 16 \ln A - 4 \gamma_E - \frac{38}{9} + \frac{8}{3} \ln(4\pi) + \frac{10}{3} \times [0.199478] \right] \right\}. \end{aligned} \qquad (3.28)$$

Here, we relabeled $\bar{\mu} \to \bar{\mu}_{UV}$ to distinguish it from other occurrences of $\bar{\mu}$, which are to be understood as factorisation scales, whereas $\bar{\mu}_{UV}$ is a genuine renormalisation scale. When merging the perturbative and lattice data, we will fix $g^2$ to the values that correspond to those chosen in the various lattice ensembles, as explained in appendix D.

## 3.4 $Z_{g,L-T}^{(c)}$ contribution

In this subsection, differently from the previous one, we rely on real-time methods: as detailed in appendix C.2.2, standard Euclidean techniques are not applicable for this class of integrals, due to the $(q^- - i\epsilon)^2$ denominator in eq. (3.20). We thus set up everything so that the integration over $K$, the loop momenta running through the self-energy blob is done in four dimensions. We choose this scheme to avoid having to consider the $D$-dimensional equivalents of the one-loop self-energies listed in appendix B. The integration over $Q$, the momentum picked up during the forward scattering, is instead done in $D = 4 - 2\epsilon$ dimensions. As we shall see, this contribution will present divergences showing sensitivity to hard, collinear and soft regions of phase space. When dealing with these regions one then needs to adopt the same scheme, as we shall do in section 4 for the $L-T$ part of the soft contribution.

<footnote n="2">²We are following the notation originally introduced in [7]; arguably, moving $g^2$ into the definition of $Z_g$ would be preferable, as it would make this operator scheme-independent.</footnote>

<footnote n="3">³The fermionic piece of the expression below has also been determined using real-time methods such as those employed in the next subsection, see [25]. The two determinations agree within numerical uncertainty.</footnote>
The fermionic and bosonic $L-T$ contributions are computed in appendix C.2.1 and C.2.2, respectively. Furthermore, as the vacuum contribution of the longitudinal and transverse self-energies is identical by Lorentz invariance, no charge renormalisation is needed in this case; indeed, the counterterm (3.27) was fully absorbed by the $(q^-)^2$ contribution.

Our result then comes from the sum of eqs. (C.69) and (C.78), viz.

$$ Z_{g,L-T,n \neq 0}^{(c)} = \frac{m_D^2}{(4\pi)^2} \left\{ \begin{aligned} & \frac{2}{\epsilon^2} + \frac{2}{\epsilon} \left( \ln \frac{\bar{\mu}e^{\gamma_E - \frac{1}{2}}}{2\pi T} + \ln \frac{\bar{\mu}A^{12}e^{-\frac{1}{2}}}{8\pi T} \right) \\ & + 2 \ln^2 \frac{\bar{\mu}e^{\gamma_E - \frac{1}{2}}}{2\pi T} - 4\gamma_1 + \frac{\pi^2}{4} - 2\gamma_E^2 + \frac{3}{2} - 2(\ln 2 - 1)^2 \\ & + 2 \ln^2 \frac{\bar{\mu}A^{12}e^{-\frac{1}{2}}}{8\pi T} + \frac{1}{2}\ln^2(2) - 2(\ln \zeta_2)^2 + \frac{2\zeta_2''}{\zeta_2} \\ & + \frac{1}{2}\mathrm{Li}_2\left(\frac{2}{3}\right) + \frac{1}{4}\mathrm{Li}_2\left(\frac{3}{4}\right) - \frac{1}{2} - \frac{1}{4}\ln(3)\ln\frac{4}{3} \end{aligned} \right\} \\ & + \frac{g^2 C_A T^2}{48\pi^2} \left\{ 0.145277 - 5.13832 \right\} \\ & - \frac{g^2 T_F N_F T^2}{48\pi^2} \left\{ \frac{2 \ln 2}{\epsilon} + 4 \ln(2) \ln \frac{\bar{\mu}A^{12}e^{-\frac{1}{2}}}{8\pi T} + 2 \ln^2(2) + 0.438129 - 1.97382 \right\}, \quad (3.29) $$

where $\zeta_s = \zeta(s)$ is the Riemann zeta function, $(\ln \zeta_s)' = \zeta'(s)/\zeta(s)$ and $\gamma_1$ is the first Stieltjes constant. For the $C_A$-proportional gauge contribution, we have again subtracted off the zero-mode associated with the $Q$ integration. Its value is given in eq. (4.2). The parts that have been obtained numerically have not been summed, to highlight their origin from the various contributions listed in appendix C.2.1 and C.2.2. What remains contains a $1/\epsilon^2$ pole, which signals the presence of double logarithmic divergences not of vacuum origin. These divergences originate from a *pole term* and a *cut term*. The former denotes the contribution to eq. (3.20) where the discontinuity implied in the “retarded minus advanced” spectral structure comes from the pole in the $1/Q^2$ propagator. The latter instead acquires spectral weight from the imaginary parts of the polarisation tensor.

To clarify the physical origin of these divergences, we now inspect these divergent terms, deferring their precise evaluation to appendix C.2.1 and C.2.2. For what concerns the pole term, it arises from eq. (C.53) and its bosonic equivalent, yielding

$$ Z_{g,L-T,\text{pole}}^{(c)} = -\frac{m_D^2}{2} \int_Q \frac{1+n_B(|q^0|)q_\perp^2}{(q^- - i\epsilon)^2 q^2} 2\pi\delta(Q^2) = -\frac{m_D^2}{2} \int_q \frac{1+2n_B(q)q_\perp^2+2q_z^2}{2q^3 q_\perp^2}. \quad (3.30) $$

While the vacuum term would be scale-free and vanishing in DR, it physically corresponds to a UV divergence for $q^z \to \infty$, which is accompanied by a *collinear divergence*. This shows up for $q^z \gg q_\perp \to 0$, i.e.<fnref n="4" />

$$ Z_{g,L-T,\text{pole div}}^{(c)} = -m_D^2 \int_q \frac{1}{2|q^z|} \frac{1}{q_\perp^2} \Big|_{q^z \gg q_\perp} . \quad (3.31) $$

While we leave the proper analysis of this divergence to future work, we think it is safe to speculate that it will be absorbed by properly dealing with the small $q_\perp$ and large $q^z$ regions.

<footnote n="4">⁴There is in addition a power-law divergent zero-mode term from the $T/q^0$ infrared limit of the Bose-Einstein distribution. This is not unexpected: it corresponds to the $O(m_D^2)$ term in the $q \gg m_D$ expansion of eq. (2.12). It too needs subtracting, as it is already included there. DR takes care of this automatically.</footnote>
The former will require systematic resummation of the asymptotic mass and of multiple soft scattering effects, likely relating this double-logarithmic divergence to the ones observed in the transverse-momentum broadening coefficient in [30, 31] and systematised for a weakly-coupled QGP in [32].<fnref n="5" /> The large-$q^z$ region is instead likely related to the emergence of the $q^z \sim E$ region, with $E \gg T$ the large light-cone energy of the hard jet parton, thus relating this divergence to the possible loop corrections to the matching coefficients in eq. (1.1).

Let us now inspect the cut term. As detailed in appendix C.2.1 and C.2.2, the divergence arises again from the vacuum part, i.e. $\theta(q^0)$, in the $1 + n_B(q^0)$ prefactor. As explained there, it is convenient to make the integrand explicitly even in $q^0$, leading to

$$Z_{g,L-T,\text{cut div}}^{(c)} = \frac{g^2}{16\pi} \int_0^\infty dk f(k) \int_0^\infty \frac{dq^0}{2\pi} \int_{\mathbf{q}} \frac{q_\perp^2}{q^5} \left[ \frac{1}{(q^0 - q^z - i\epsilon)^2} + \frac{1}{(q^0 + q^z + i\epsilon)^2} \right] \\ \times \theta(q - 2k) \left[ (3(2k - q^0)^2 - q^2) \theta(q + 2k - q^0)\theta(q^0 - q) \right. \\ \left. - (3(2k + q^0)^2 - q^2) \theta(q^0 - q + 2k)\theta(q - q^0) \right], \quad (3.32)$$

where

$$f(k) = 2T_F N_F n_F(k) + C_A n_B(k), \quad (3.33)$$

and $n_F(k) = (\exp(k/T) + 1)^{-1}$ is the Fermi-Dirac distribution. Eq. (3.32) arises from appropriately combining eqs. (C.55) and (C.65), together with their bosonic counterparts like eq. (C.74). We only consider here the $\theta(q - 2k)$ part of eqs. (C.65) and (C.74) as it is the divergent one.

Let us consider for illustration the simplest part of the above, the $q^2$-proportional piece in $3(2k - q^0)^2 - q^2)$ in the time-like slice $q^0 > q$. The main message will not be altered by including the other terms. It reads

$$\int_0^\infty dk f(k) \int_{\mathbf{q}} \int_q^{q+2k} \frac{dq^0 q_\perp^2}{2\pi q^3} \left[ \frac{1}{(q^0 - q^z)^2} + \frac{1}{(q^0 + q^z)^2} \right] \theta(q - 2k) \\ = 2 \int_0^\infty dk k f(k) \int_{\mathbf{q}} \frac{2(kq + q_z^2) + q_\perp^2}{\pi q^3 (q_\perp^2 + 4kq + 4k^2)} \theta(q - 2k) \quad (3.34)$$

$$= \int_0^\infty dk k f(k) \int_{q_\perp} \frac{4q_\perp^2}{\pi^2(q_\perp^2 + 4k^2)^2} \left[ \ln \frac{q_\perp}{2k} + \mathcal{O}(1) + \mathcal{O}(k/q_\perp) \right], \quad (3.35)$$

where the arising logarithmic term is what causes the double pole in DR. Its $k$-dependent argument is what is responsible for the different arguments of the logarithms in the fermionic and bosonic case, see the $N_F$-proportional $\ln 2$ term in eq. (3.29). Moreover, the $\ln(q_\perp/(2k))$ must come from two regions of the $q^z$ integral in eq. (3.34), once we take the asymptotic UV expansion $q_\perp \gg k$. Indeed, we have a first region for $q^z \sim q_\perp \gg k$ and a second region for $q_z \sim q_\perp^2/k \gg q_\perp$ — we are exploiting evenness in $q^z$ to consider twice the positive range only. We then have, by appropriately expanding the integrand in the two regions

<footnote n="5">⁵It also remains to be understood how such a long-duration collinear process could be used as an ingredient in the determination of an equally long-duration (formation time) medium-induced radiation. This question too is left for future work.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** EQCD equivalent of diagram (c). The labelling is identical to that from figure 1 except that here, curly lines represent spatial gauge bosons, double lines are adjoint Wilson lines and solid lines are adjoint scalars.</figcaption></figure>

and separating them by a cutoff $\mu$

$$ \int_0^\mu \frac{dq^z}{\pi} \frac{2q_z^2+q_\perp^2}{\pi q^3 q_\perp^2} = \frac{\ln \frac{4\mu^2}{q_\perp^2} - 1}{\pi^2 q_\perp^2}, \quad \text{and} \quad \int_\mu^\infty \frac{dq^z}{\pi} \frac{2}{\pi q^z (q_\perp^2 + 4kq^z)} = \frac{\ln \frac{q^4}{16k^2\mu^2}}{\pi^2 q_\perp^2}, \qquad (3.36) $$

which reproduces the $q_\perp \gg k$ limit of eq. (3.35).

This then clarifies that the cut term has a double-logarithmic UV divergence, which we can again speculate to be related to the emergence of the $q^+ \sim E$ region. Its proper handling is also left to future work. We conclude by remarking that in the real-time calculation of $Z_{g,q^{-2}}^{(c)}$ (cf. [25]) there is a cancellation between a double-log akin to that in eq. (3.35) and an opposite one arising from a pole contribution which is very different from the simple one in eq. (3.31). The collinear sensitivity of the latter and its absence in $Z_{g,q^{-2}}^{(c)}$ — see [25] — might explain why no such cancellation happens for $Z_{g,L-T}^{(c)}$.

# 4 Cancellation of IR QCD divergences with UV EQCD divergences

As we explained around eq. (2.18), the EQCD evaluation in [6] relied on a subtraction of the $c_2$-proportional UV-log divergent term from NLO perturbative EQCD in eq. (2.18), which introduced a logarithmic dependence on the $L_{min}$ regulator. The Q zero modes in the perturbative QCD evaluation of the previous section, on the other hand, introduced an IR log-divergence, viz.

$$ Z_{g,q^{-2},B,n=0}^{(c)} = \frac{3g^2 N_c T^2}{64\pi^2} \left[ -\frac{1}{\epsilon} - 4 \ln \left( \frac{\bar{\mu}}{2T} \right) - \frac{5}{3} \right], \qquad (4.1) $$

$$ Z_{g,L-T,B,n=0}^{(c)} = \frac{g^2 C_A T^2}{32\pi^2} \left[ -\frac{1}{\epsilon} - 2 \ln \left( \frac{\bar{\mu}}{4T} \right) + 2 \right], \qquad (4.2) $$

which have been obtained from eqs. (C.17) and (C.71).

Our goal in this section is to add back a DR version of the NLO subtraction in eq. (2.18). The resulting UV poles in $1/\epsilon$ should cancel against those in eqs. (4.1) and (4.2). To this end, we note that [6] used Feynman gauge for the 3D gauge field of EQCD. From the arguments of section 3.1 we then know that this $c_2$-proportional UV divergence must come from the EQCD analogue of diagram (c), shown in figure 2. This suggests that

$$ Z_g^{\text{match}} = Z_{g,q^{-2},B,n=0}^{(c)} + Z_{g,q^{-2},c_2}^{3d(c)} + Z_{g,L-T,B,n=0}^{(c)} + Z_{g,L-T,c_2}^{3d(c)}, \qquad (4.3) $$

where $Z_{g,q^{-2},c_2}^{3d(c)}$ is the $(q^-)^2$ and $Z_{g,L-T,c_2}^{3d(c)}$ the $L-T$ piece of the UV divergence, evaluated in DR. When adding eq. (4.3) to (2.18), we will undo the $c_2$ subtraction and thus arrive at

$$ Z_g^{\text{non-pert class}} = Z_g^{3d} |^{\text{merge}} + Z_g^{\text{match}}, \qquad (4.4) $$

which is a finite, $L_{min}$-independent result for the non-perturbative evaluation of the classical
contribution. The best-available determination of $Z_g$ is then

$$Z_g^{\text{non-pert class+NLO T}} = Z_g^{\text{LO}} + \delta Z_g^{\text{non-pert class+NLO T}}, \quad (4.5)$$

$$\delta Z_g^{\text{non-pert class+NLO T}} = Z_g^{\text{non-pert class}} + Z_{g,q^{-2}}^{(c),\text{IR/UV safe}} + Z_{g,L-T,n \neq 0}^{(c)}, \quad (4.6)$$

where $\delta Z_g^{\text{non-pert class+NLO T}}$ incorporates the classical modes to all orders and the thermal modes to NLO. The final two terms in eq. (4.6) have been given in eqs. (3.28) and (3.29).

Let us now begin the determination of $Z_{g,q^{-2},c_2}^{3d(c)}$ and $Z_{g,L-T,c_2}^{3d(c)}$, starting from the former. With respect to the decomposition defined in eq. (2.16), we first take the $\langle BB \rangle$ contribution from appendix A.3 of [6] (cf. eq. (A.14))

$$\langle BB \rangle_{\text{div}}^{(c)} = \frac{1}{4} \int_{\mathbf{q}} \frac{e^{iLq_z}}{q^4} (q_{\perp}^2 + (d-1)q_z^2) \delta^{ab} \left( - \Pi_{A_i^a A_j^b \text{div}}^T(q) \right), \quad (4.7)$$

where $\Pi_{A_i^a A_j^b \text{div}}^T$ is the EQCD transverse gluon polarisation tensor in the massless limit. As we are investigating a UV divergence on the EQCD side, the $q \gg m_D$ limit is the appropriate one. It can be extracted either from the $m_D \to 0$ limit of eq. (A.10) in [6] or from the zero-mode contribution to eq. (B.4). We find

$$\Pi_{A_i^a A_j^b \text{div}}^T(q) = \delta^{ab} \frac{3g_{3d}^2 C_A}{(4\pi)^{\frac{d}{2}}} \frac{\Gamma(1-\frac{d}{2})\Gamma(\frac{d}{2})^2}{\Gamma(d-1)} \mu^{3-d} q^{d-2}, \quad (4.8)$$

$$\begin{align}
& d=3-2\epsilon \nonumber \\
& \equiv -\delta^{ab} \frac{3g_{3d}^2 C_A q}{16}, \tag{4.9}
\end{align}$$

where we have denoted the $d=3$ limit for future use.

The second term in eq. (4.7), proportional to $(d-1)q_z^2$, corresponds to the $Z_{g,q^{-2},c_2}^{3d(c)}$ contribution. To obtain it, we must then integrate in $L$ up to the cutoff while performing the $q$ integration in DR, so that its $Z_g^{\text{match}}$ contribution undoes the corresponding piece of the $c_2$ subtraction. Looking at eq. (2.16) and using eq. (4.8), we then find

$$\begin{align}
Z_{g,q^{-2},c_2}^{3d(c)} &= -\frac{4T}{d_A} \int_0^{L_{\min}} dL \, L \, \langle BB \rangle_{q^{-2} \text{div}} = \frac{3g_{3d}^2 C_A T}{128\pi^d} \frac{(d-1)(2d-5)}{(d-3)} \frac{\Gamma(\frac{d}{2}-2)\Gamma(\frac{d}{2}-1)}{(L_{\min}\mu)^{2d-6}} \nonumber \\
&= \frac{3g_{3d}^2 C_A T}{64\pi^2} \left[ \frac{1}{\epsilon} - 7 + 4\gamma_E + 4\ln(L_{\min}\bar{\mu}) \right]. \tag{4.10}
\end{align}$$

Since $g_{3d}^2 = g^2 T + O(g^4)$, the $1/\epsilon$ term cancels against the one from eq. (4.1).

As for the $L-T$ piece, we need the $q_{\perp}^2$-proportional piece of eq. (4.7) together with the contribution from the $\langle EE \rangle$ correlator. We extract this from eq. (A.8) of [6] and expand towards the UV to arrive at

$$\begin{align}
\Pi_{\Phi^a \Phi^b}(q) &= \delta^{ab} \frac{g_{3d}^2 C_A m_D}{16\pi} \left( -4 - \frac{6C_F - C_A}{6C_A} \frac{\lambda_E}{g_{3d}^2} - 8 \frac{q^2 - m_D^2}{q m_D} \arctan \frac{q}{m_D} \right) \nonumber \\
&\quad \approx -\delta^{ab} \frac{g_{3d}^2 C_A q}{4} \equiv \Pi_{\Phi^a \Phi^b \text{div}}(q). \tag{4.11}
\end{align}$$

As the corresponding IR pole has been obtained using the mixed scheme described at the beginning of section 3.4, we stick to the same scheme here and have evaluated eq. (4.11) in $d=3$. From eq. (A.13) of [6], we then write

$$\langle EE \rangle_{\text{div}}^{(c)} = \frac{1}{4} \int_{\mathbf{q}} \frac{e^{iLq_z}}{q^4} q_{\perp}^2 \delta^{ab} \left( - \Pi_{\Phi^a \Phi^b \text{div}}^T(q) \right). \quad (4.12)$$
To get the full $\Pi_L - \Pi_T$ EQCD contribution, we need to add to this its $\langle BB \rangle$ counterpart. It arises from the $q_\perp^2$-proportional piece in eq. (4.7) with the $d=3$ polarisation tensor in eq. (4.9). We then evaluate the Fourier transform in DR and carry out the $L$ integration, as in eq. (4.10). This yields — see again eq. (2.16)

$$ 
\begin{aligned}
Z_{g \ L-T \ c_2}^{3d \ (c)} &= \frac{4T}{d_A} \int_0^{L_{\min}} dL \, L \left( \langle EE \rangle_{\text{div}}^{(c)} - \langle BB \rangle_{\text{div}}^{(c)} \right) = \frac{g_{3d}^2 C_A T}{(4\pi)^{\frac{d+1}{2}}} \frac{\Gamma(\frac{d+1}{2})}{(3-d)} \left( \frac{L_{\min} \mu}{2} \right)^{3-d} \\
&= \frac{g_{3d}^2 C_A T}{32\pi^2} \left[ \frac{1}{\epsilon} - 1 + 2\gamma_E + 2 \ln \frac{L_{\min} \bar{\mu}}{2} \right],
\end{aligned}
\quad (4.13) $$

which, reassuringly contains a $1/\epsilon$ pole that cancels exactly against the one from eq. (4.2). We note that, if the Fourier transforms in eqs. (4.7) and (4.12) are carried out in three dimensions, the former vanishes, whereas the latter gives $\langle EE \rangle_{\text{div } d=3}^{(c)} = C_A d_A g_{3d}^2 / (4\pi L)^2$, from which the value of $c_2$ in eq. (2.19) is extracted. We have split our calculation in a $(q^-)^2$ part and a $L-T$ part, rather than a $\langle EE \rangle$ and a $\langle BB \rangle$ one, to check the explicit cancellation of $1/\epsilon$ poles in each of the former two parts and to account for our mixed scheme in the $L-T$ contribution.

We can finally insert eqs. (4.1), (4.2), (4.10) and (4.13) into eq. (4.3), together with the LO EQCD matching $g_{3d}^2 = g^2 T$, to find

$$ Z_{g}^{\text{match}} = \frac{g^2 C_A T^2}{8\pi^2} \left[ 2 \ln(2L_{\min}T) + 2\gamma_E - 3 \right]. \quad (4.14) $$

Note that in computing $Z_g^{\text{non-pert class}}$, the $\ln L_{\min}$ above cancels against the $-\int \frac{dL}{L} c_2 g_{3d}^2$ term in eq. (2.18), as we originally set to show.

We conclude by noting that eq. (4.14) is somewhat arbitrary: a firm necessity of the cancellation of the $L_{\min}$ dependence of $Z_g^{3d}|_{\text{merge}}$ through the IR divergence of the thermal scale contribution is the $\ln(L_{\min}T)$ term. The non-logarithmic parts, on the other hand, can freely be shuffled between eq. (4.14) and eqs. (3.28) and (3.29) by redefining the scheme used to separate the IR-divergent part from the remainder in section 3. This is of course irrelevant once these three equations are summed in eq. (4.6). However, in the next section we shall not include the $Z_{g \ L-T \ n \neq 0}^{(c)}$ contribution, as it contains extra divergences signaling sensitivities to new regions of phase space we have left to future work. Hence, to present a first estimate of the impact of neglecting $Z_{g \ L-T \ n \neq 0}^{(c)}$, let us determine what would happen to eq. (4.14) if we were to omit the finite parts associated with the $L-T$ contribution. By dropping them from eqs. (4.2) and (4.13) we find

$$ Z_{g, \text{no } L-T \text{ finite parts}}^{\text{match}} = Z_{g}^{\text{match}} - \frac{g^2 C_A T^2}{32\pi^2} \left[ 1 + 2 \ln 2 + 2\gamma_E \right]. \quad (4.15) $$

The numerical impact of eqs. (4.14) and (4.15) will be discussed in the next section.

# 5 Summary and discussion

In section 3, we determined the contribution to $Z_g$ of the temperature scale at order $g^2$: it is given by eqs. (3.28) and (3.29) for the two contributions in which it has naturally been split; see eq. (3.18). In section 4, we analyzed the IR divergences of the $T$-scale contribution and showed how they cancel against the UV divergences from EQCD found in [6] when they
are consistently evaluated in the same DR scheme. In eq. (4.4), we explain how to add the lattice result of [6] with the finite remainder, eq. (4.14), of the cancellation of the divergences, lifting the artificial dependence on the UV cutoff introduced in [6].

To evaluate the numerical impact of our work, we should in principle use our best available result, summing the NLO T-scale contribution with the all-order non-perturbative evaluation. This is given by eq. (4.6), which, after inserting all ingredients, gives rise to

$$
\begin{equation}
\begin{aligned}
\delta Z_g^{\text{non-pert class+NLO T}} = & Z_g^{3d|\text{merge}} + \frac{g^2 C_A T^2}{8\pi^2} \left[ 2\ln(2L_{\text{min}T}) + 2\gamma_E - 3 \right] \\
& + \frac{g^2 T^2}{96\pi^2} \left\{ \left( \frac{11}{3}C_A - \frac{4}{3}T_F N_F \right) 2\ln \frac{\bar{\mu}_{\text{UV}}}{T} \right. \\
& \left. \qquad + C_A \left[ -80\ln A - 4\gamma_E + \frac{94}{9} + \frac{32}{3}\ln\pi - \frac{8}{3}\ln 2 + \frac{4}{3} \times [0.299378] \right] \right. \\
& \left. \qquad + T_F N_F \left[ 16\ln A - 4\gamma_E - \frac{38}{9} + \frac{8}{3}\ln(4\pi) + \frac{10}{3} \times [0.199478] \right] \right\} \\
& + \frac{m_D^2}{(4\pi)^2} \left\{ \frac{2}{\epsilon^2} + \frac{2}{\epsilon} \left( \ln \frac{\bar{\mu}e^{\gamma_E - \frac{1}{2}}}{2\pi T} + \ln \frac{\bar{\mu}A^{12}e^{-\frac{1}{2}}}{8\pi T} \right) \right. \\
& \left. \qquad + 2\ln^2 \frac{\bar{\mu}e^{\gamma_E - \frac{1}{2}}}{2\pi T} - 4\gamma_E + \frac{\pi^2}{4} - 2\gamma_E^2 - 1 + 4\ln 2 \right. \\
& \left. \qquad + 2\ln^2 \frac{\bar{\mu}A^{12}e^{-\frac{1}{2}}}{8\pi T} - \frac{3}{2}\ln^2(2) - 2(\ln \zeta_2)^{\prime 2} + \frac{2\zeta_2^{\prime\prime}}{\zeta_2} \right. \\
& \left. \qquad + \frac{1}{2}\text{Li}_2\left(\frac{2}{3}\right) + \frac{1}{4}\text{Li}_2\left(\frac{3}{4}\right) - \frac{1}{4}\ln(3)\ln\left(\frac{4}{3}\right) \right\} \\
& + \frac{g^2 T^2}{48\pi^2} \left\{ C_A \left[ 0.145277 - 5.13832 \right] \right. \\
& \left. \qquad - T_F N_F \left[ \frac{2\ln 2}{\epsilon} + 4\ln(2)\ln \frac{\bar{\mu}A^{12}e^{-\frac{1}{2}}}{8\pi T} - 2\ln^2(2) + 0.438129 - 1.97382 \right] \right\}.
\end{aligned}
\tag{5.1}
\end{equation}
$$

The one outstanding roadblock is that our evaluation of the “longitudinal minus transverse” $Z_{g,L-T}^{(c)} \neq 0$ contribution in eq. (3.29) contains yet uncancelled double-logarithmic divergences from the sensitivity to the collinear ($q_\perp \sim m_\infty$) and hard ($q^+ \sim E \gg T$) regions of phase space. We have left their understanding and cancellation to future work. Here, we present a numerical analysis of the impact of the other corrections, namely the $(q^-)^2$ contribution, $Z_{g,q^{-2}}^{(c),\text{IR/UV safe}}$, given in eq. (3.28) and the classical contribution determined from lattice and perturbative EQCD in [6] complemented by the cancellation of the logarithmic sensitivity achieved through eqs. (4.3) and (4.14). More precisely, let us define

$$
\begin{align}
\delta Z_g^{\text{finite}} &\equiv Z_g^{\text{non-pert class}} + Z_{g,q^{-2}}^{(c),\text{IR/UV safe}} \tag{5.2} \\
&= Z_g^{3d|\text{merge}} + \frac{g^2 C_A T^2}{8\pi^2} \left[ 2 \ln \frac{T}{2g_{3d}^2} + 2\gamma_E - 3 \right] \nonumber \\
&\quad + \frac{g^2 T^2}{96\pi^2} \left\{ \left( \frac{11}{3} C_A - \frac{4}{3} T_F N_F \right) 2 \ln \left( 4\pi e^{-\gamma_E} \exp \frac{-N_c + 4 \ln(4) N_F}{22 N_c - 4 N_F} \right) \right. \nonumber \\
&\quad \left. + C_A \left[ -80 \ln A - 4\gamma_E + \frac{94}{9} + \frac{32}{3} \ln \pi - \frac{8}{3} \ln 2 + \frac{4}{3} \times [0.299378] \right] \right. \nonumber \\
&\quad \left. + T_F N_F \left[ 16 \ln A - 4\gamma_E - \frac{38}{9} + \frac{8}{3} \ln(4\pi) + \frac{10}{3} \times [0.199478] \right] \right\}, \tag{5.3}
\end{align}
$$
<table><caption>**Table 1.** Values for our main results. $\delta Z_{g,Eq. (4.15)}^{finite}$ is defined by replacing eq. (4.14) with eq. (4.15) in $Z_{g}^{non-pert class}$ and $Z_{g,LO}^{3d}$ is given by eq. (2.12). See appendix D for details on the renormalisation scale and the values of the coupling. We recall that these values are to be added to the LO value $Z_{g}^{LO} = T^2/6$, as per eq. (4.5).</caption><thead><tr><th>T</th><th>Z<sub>g</sub><sup>non-pert class</sup><br>T<sup>2</sup></th><th>δZ<sub>g</sub><sup>finite</sup><br>T<sup>2</sup></th><th>δZ<sub>g,Eq. (4.15)</sub><br>T<sup>2</sup></th><th>Z<sub>g,LO</sub><sup>3d</sup><br>T<sup>2</sup></th></tr></thead><tbody><tr><td>250 MeV</td><td>-0.513(138)(45)(7)</td><td>-0.340(138)(45)(7)</td><td>-0.465(138)(45)(7)</td><td>-0.376</td></tr><tr><td>500 MeV</td><td>-0.619(99)(39)(3)</td><td>-0.491(99)(39)(3)</td><td>-0.584(99)(39)(3)</td><td>-0.324</td></tr><tr><td>1 GeV</td><td>-0.462(71)(9)(7)</td><td>-0.356(71)(9)(7)</td><td>-0.430(71)(9)(7)</td><td>-0.305</td></tr><tr><td>100 GeV</td><td>-0.327(16)(5)(2)</td><td>-0.274(16)(5)(2)</td><td>-0.310(16)(5)(2)</td><td>-0.223</td></tr></tbody></table>

which is obtained by summing eqs. (3.28) and (4.3). We also set $\bar{\mu}_{UV}$ to eq. (D.1), in keeping with the renormalisation scale used in the EQCD matching, and $L_{min} = 1/(4g_{3d}^2)$, the value used in $Z_{g}^{3d}|^{merge}$ in [6]. To further motivate our choice of looking, for the time being, at this partial subset of $O(g^2)$ corrections, let us remark that $Z_{g,q^{-2}}^{(c) IR/UV safe}$ comes from the same $(q^{-})^2$ structure as the LO term — compare eq. (2.3) and (3.19). It is thus fundamentally different from $Z_{g,L-T,n \neq 0}^{(c)}$, which indeed contains a new structure and new divergences.

In table 1, we collect the numerical values for our main results. We refer to appendix D for the coupling prescriptions and the EQCD data for $Z_{g}^{3d}|^{merge}$ from [6]. $\delta Z_{g,Eq. (4.15)}^{finite}$ is defined by replacing eq. (4.14) with eq. (4.15) in $Z_{g}^{non-pert class}$, thus giving a first, very rough estimate of the error associated with the lack of inclusion of the $L-T$ contribution. The quoted uncertainties are inherited from [6]. They come from statistical errors in lattice EQCD (first bracket), from the integration of the IR tail in eq. (4.3) (second bracket) and from the quadrature of the lattice data (third bracket).

As we see, our results for $\delta Z_{g}^{finite}$ are negative and, in magnitude, larger than $Z_{g}^{LO} = T^2/6$ (see eq. (2.3)), to which they are to be added, as per eq. (4.5). With the exception of the highest temperature, within errors they are however compatible with the $O(g)$ correction only, i.e. the leading-order contribution in perturbative EQCD given by $Z_{g,LO}^{3d}$ in the last column of the table. Recall that $Z_{g,LO}^{3d}$ is obtained in eq. (2.12), which in turn clarifies why a negative contribution arises: it comes from removing the incorrect, unscreened IR of bare QCD and replacing it with the correct, screened IR of EQCD. As the screened contribution is necessarily smaller than the bare unscreened one, negativity sets in.

This large negative contribution, resulting in an overall negative $Z_g$, may appear problematic. However, until the divergences of the $L-T$ contributions are addressed by properly including collinear and hard modes, we refrain from making definite statements on the convergence of EQCD factorisation for this observable. Already by simply reshuffling the finite parts associated with the cancellation of IR and UV divergences in the $L-T$ channel — the $\delta Z_{g,Eq. (4.15)}^{finite}$ column — we see a significant variation. Fully understanding these extra regions of phase space may give rise to even larger contributions, due to potentially large double logarithms.

These potentially enhanced double logarithms would relate the $T$-scale physics in eqs. (5.1) and eq. (3.29) with hard, $q^+ \sim E \gg T$ physics related to the EFT matching implicit in eq. (1.1), as we commented there. They would also include collinear physics, with effects such as the quantum-mechanical interference of multiple soft scatterings, as we mentioned
in section 3.4 after eq. (3.31). At the moment we do not speculate further on these matters and leave them to future work.

As a final remark, we comment on the expansion undergirding the factorisation of the asymptotic mass into the operators $Z_g$ and $Z_f$ described by eq. (1.1). It requires the parton energy $E$ to be $E \gg T$ when the $T$-scale contribution to these operators is considered. If the soft contribution is considered, the expansion becomes $E \gg m_D$, so that, in the weak-coupling regime, it remains applicable for $E \sim T$. This also holds for eqs. (4.4) and (4.14), as the latter is obtained from the soft limit of the $T$-scale contribution.<fnref n="6" /> This has then the very important consequence that our values for $Z_g^{\text{non-pert class}}$, tabulated in table 1, are valid also when considering partons with thermal energies. This then makes them applicable for corrections to observables such as the thermal photon rate or transport coefficients, where they could be used to complement the existing NLO determinations [33–35].

## Acknowledgments

We thank Guy Moore for useful conversations. JG, PS, NS and EW acknowledge support by a PULSAR grant from the Région Pays de la Loire. JG acknowledges support by the Agence Nationale de la Recherche under grant ANR-22-CE31-0018 (AUTOTHERM). PS acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) through the CRC-TR 211 ‘Strong-interaction matter under extreme conditions’ — project number 315477589 — TRR 211.

## A Conventions

Our sign for the covariant derivative is

$$D_{\mu} = \partial_{\mu} - igA_{\mu},$$

which fixes the sign of the three gluon vertex to be positive. Moreover, we employ the “mostly minus” $(+, -, -, -)$ metric. Uppercase letters denote four-momenta, lowercase letters the modulus of the three-momenta. When using light-cone coordinates

$$p^{+} = \frac{p^{0} + p^{z}}{2} = \bar{v} \cdot P, \quad p^{-} = p^{0} - p^{z} = v \cdot P, \quad P \cdot Q = p^{+}q^{-} + p^{-}q^{+} - p_{\perp} \cdot q_{\perp}, \quad (A.1)$$

where the two light-like reference vectors are defined as

$$\bar{v}^{\mu} \equiv \frac{1}{2}(1, 0, 0, -1), \quad v^{\mu} \equiv (1, 0, 0, 1).$$

This asymmetric convention for the $+$ and $-$ components of the light-cone coordinates has two advantages: it has unitary Jacobian, i.e. $dp^0dp^z = dp^+dp^-$, in scalings where $p^- \ll p^+$ it then implies $p^0 \approx p^z \approx p^+$.

<footnote n="6">⁶We remark that, even in this regime, our results remain different from those of [8–11]. These papers determine the two-loop and power-suppressed one-loop contributions from the hard thermal scale to the self-energy for soft external modes, to then take the asymptotic mass limit. It then corresponds to first expanding for $P \ll Q_i \sim T$, where $P$ is the external momentum and $Q_i$ the loop momenta, to then take the $P \gg m_D$ limit. Our calculation, on the other hand, corresponds to taking a $P \gg Q_i \sim T$ expansion to then take the $P \to T \gg m_D$ limit from above. Indeed, our calculation in this limit sees IR divergences given by eqs. (4.1)–(4.2) which are then absorbed by EQCD, whereas the results of [8–11] are finite after charge renormalization, further highlighting their difference.</footnote>
The emergent momenta integrals are often regularised in $d = D - 1 = 3 - 2\epsilon$, employing the $\overline{\text{MS}}$ renormalisation scale $\bar{\mu}^2 = 4\pi e^{-\gamma_E} \mu^2$. Denoting spatial momenta as $\mathbf{p} = (\mathbf{p}_{\perp}, p_z)$, this allows to use the shorthand

$$ \begin{aligned} \int_{P_E} \equiv T \sum_{n \in \mathbb{Z}} \int_{\mathbf{p}} \Big|_{p_0=2n\pi T}, & \qquad \int_P \equiv \mu^{3-d} \int \frac{d^D P}{(2\pi)^D}, & \qquad \int_{\mathbf{p}} \equiv \mu^{3-d} \int \frac{d^d \mathbf{p}}{(2\pi)^d}, \\ \int_{\{P_E\}} \equiv T \sum_{n \in \mathbb{Z}} \int_{\mathbf{p}} \Big|_{p_0=(2n-1)\pi T}, & \qquad \int_{\mathbf{p}_{\perp}} \equiv \mu^{3-d} \int \frac{d^{d-1} \mathbf{p}_{\perp}}{(2\pi)^{d-1}}. & \end{aligned} \quad (A.2) $$

## B One-loop self-energies and projectors

We document in this appendix the self energies that are needed for the computation of diagram (c) in section 3. We start by providing the transverse and longitudinal projectors (with respect to the three-momentum) that are used to arrive at the decomposition eq. (3.18)

$$ P_L^{\mu\nu}(Q) = -g^{\mu\nu} + \frac{q^\mu q^\nu}{Q^2} - P_T^{\mu\nu}(Q), \quad (B.1) $$

$$ P_T^{\mu\nu}(Q) = -g^{\mu\nu} + u^\mu u^\nu - \frac{q^\mu q^\nu - (Q \cdot u)q^\mu u^\nu - (Q \cdot u)q^\nu u^\mu + (Q \cdot u)^2 u^\mu u^\nu}{q^2}, \quad (B.2) $$

$$ (u^\mu) = (1, 0, 0, 0). \quad (B.3) $$

The main ingredient for the Euclidean evaluation of the $(q^-)^2$ contribution is the Feynman-gauge transverse Euclidean polarisation tensor, which can be obtained from, e.g. [36]. It reads

$$ \begin{aligned} \Pi_T^E(Q) ={}& g^2 N_c \int_{K_E} \frac{(D-2)K_E^2 - Q_E^2 - \frac{1}{2} \frac{Q_E^2}{q^2} (Q_E^2 + 2k^2 - 2k_0^2)}{K_E^2 (Q_E - K_E)^2} \\ & - 2g^2 T_F N_F \int_{\{K_E\}} \frac{2K_E^2 - \frac{D-4}{D-2} Q_E^2 - \frac{Q_E^2}{(D-2)q^2} (Q_E^2 + 2k^2 - 2k_0^2)}{K_E^2 (Q_E - K_E)^2}. \end{aligned} \quad (B.4) $$

For the $L-T$ contribution, we need instead the difference between the retarded longitudinal and transverse part, which appears in eq. (3.20). Its gauge contribution (gluon and ghost) and its quark contribution can be read off from [37], viz.

$$ \begin{aligned} \Pi_L^R(Q)_B - \Pi_T^R(Q)_B ={}& \frac{g^2 C_A}{(4\pi)^2 q^3} \int_0^\infty dk n_B(k) \left[ -8kq(3q_0^2 - q^2) \right. \\ & \left. \qquad + Q^2 (12k^2 - q^2 + 3q_0^2) \left[ \ln \frac{q^0 + q - 2k + i\epsilon}{q^0 - q - 2k + i\epsilon} + \ln \frac{q^0 - q + 2k + i\epsilon}{q^0 + q + 2k + i\epsilon} \right] \right. \\ & \left. \qquad + 12kq^0 Q^2 \left[ 2 \ln \frac{q^0 + q + i\epsilon}{q^0 - q + i\epsilon} + \ln \frac{q^0 - q + 2k + i\epsilon}{q^0 + q + 2k + i\epsilon} - \ln \frac{q^0 + q - 2k + i\epsilon}{q^0 - q - 2k + i\epsilon} \right] \right], \end{aligned} \quad (B.5) $$

$$ \begin{aligned} \Pi_L^R(Q)_F - \Pi_T^R(Q)_F ={}& \frac{2g^2 T_F N_F}{(4\pi)^2 q^3} \int_0^\infty dk n_F(k) \left[ -8kq(3q_0^2 - q^2) \right. \\ & \left. \qquad + Q^2 (12k^2 - q^2 + 3q_0^2) \left[ \ln \frac{q^0 + q - 2k + i\epsilon}{q^0 - q - 2k + i\epsilon} + \ln \frac{q^0 - q + 2k + i\epsilon}{q^0 + q + 2k + i\epsilon} \right] \right. \\ & \left. \qquad + 12kq^0 Q^2 \left[ 2 \ln \frac{q^0 + q + i\epsilon}{q^0 - q + i\epsilon} + \ln \frac{q^0 - q + 2k + i\epsilon}{q^0 + q + 2k + i\epsilon} - \ln \frac{q^0 + q - 2k + i\epsilon}{q^0 - q - 2k + i\epsilon} \right] \right], \end{aligned} \quad (B.6) $$

where B and F stand for bosonic and fermionic.
C Diagram (c) computation: technical details

C.1 Euclidean evaluation of $Z_{g,q-2}$

We start from eq. (3.22). Upon inserting eq. (B.4) into eq. (3.22), we encounter the following class of sum-integrals

$$ \mathcal{I}_{s_0; s_1 s_2 s_3}^{\alpha_1 \alpha_2} = \oint_{Q_E K_E} \frac{1}{[q^2]^{s_0} [Q_E^2]^{s_1} [K_E^2]^{s_2} [(Q_E - K_E)^2]^{s_3}}, \quad (C.1) $$

$$ \hat{\mathcal{I}}_{s_0; s_1 s_2 s_3}^{\alpha_1 \alpha_2} = \oint_{Q_E \{K_E\}} \frac{1}{[q^2]^{s_0} [Q_E^2]^{s_1} [K_E^2]^{s_2} [(Q_E - K_E)^2]^{s_3}}, \quad (C.2) $$

where $\mathcal{I}_{0; s_1 s_2 s_3}^{00} = \mathcal{I}_{s_1 s_2 s_3}$ and fermionic momenta are denoted by curly brackets. The integrals that appear in eq. (B.4) are

$$ J_1 = \mathcal{I}_{201} = \oint_{Q_E K_E} \frac{1}{Q_E^4 K_E^2} = \mathcal{I}_1 \mathcal{I}_2, \quad (C.3) $$

$$ J_2 = \mathcal{I}_{111} = 0, \quad (C.4) $$

$$ J_3 = \mathcal{I}_{1;011} = \oint_{Q_E K_E} \frac{1}{Q_E^2 K_E^2 (\mathbf{q} - \mathbf{k})^2} = \frac{T^2}{192\pi^2} \left[ \frac{1}{\epsilon} - 4 \ln \frac{32\pi T}{\mu A^{12}} + 6 \right], \quad (C.5) $$

$$ J_4 = \mathcal{I}_{1;101} = \oint_{Q_E} \frac{1}{q_0^2} \left[ \frac{1}{q^2} - \frac{1}{Q_E^2} \right] \oint_{K_E} \frac{1}{K_E^2} = \frac{2}{(d-2)} \mathcal{I}_1 \mathcal{I}_2, \quad (C.6) $$

$$ J_5 = \mathcal{I}_{1;111}^{01} = \frac{T^2}{288\pi^2} \left[ \frac{1}{\epsilon} + 4 \ln \frac{\mu A^3}{4\pi T} + 3\gamma_E + \frac{95}{12} - 0.299378 \right], \quad (C.7) $$

and the fermionic counterparts $\hat{J}_i$ obtained by summing over odd $\{K_E\}$. We used the one-loop bosonic and fermionic master integrals

$$ \mathcal{I}_s = \oint_{P_E} \frac{1}{[P_E^2]^s} = \mu^{3-d} 2T \frac{[2\pi T]^{d-2s}}{(4\pi)^{\frac{d}{2}}} \frac{\Gamma(s-\frac{d}{2})}{\Gamma(s)} \zeta_{2s-d}, \quad (C.8) $$

$$ \hat{\mathcal{I}}_s = \oint_{\{P_E\}} \frac{1}{[P_E^2]^s} = (2^{2s-d} - 1) \mathcal{I}_s, \quad (C.9) $$

where $\zeta_s = \zeta(s)$ is the Riemann zeta function. $J_2$ vanishes identically, as shown with integration-by-parts (IBP) methods [38]. The integrals $J_3$ in eq. (C.5) and $J_5$ in eq. (C.7) need to be evaluated independently due to absence of collinearity in the integral [39]. Details on their evaluation are provided below in appendix C.1.1.

The fermionic counterparts read

$$ \hat{J}_1 = \hat{\mathcal{I}}_{201} = \oint_{Q_E \{K_E\}} \frac{1}{Q_E^4 K_E^2} = \hat{\mathcal{I}}_1 \mathcal{I}_2, \quad (C.10) $$

$$ \hat{J}_2 = \hat{\mathcal{I}}_{111} = 0, \quad (C.11) $$

$$ \hat{J}_3 = \hat{\mathcal{I}}_{1;011} = \oint_{\{Q_E K_E\}} \frac{1}{Q_E^2 K_E^2 (\mathbf{q} - \mathbf{k})^2} = \frac{-T^2}{96\pi^2} \left[ \frac{1}{\epsilon} - 4 \ln \frac{8\pi T}{\mu A^{12}} + \ln 4 \right], \quad (C.12) $$
$$ \hat{J}_4 = \hat{\mathcal{I}}_{1;101} = \sum_{Q_E} \int \frac{1}{q_0^2} \left[ \frac{1}{q^2} - \frac{1}{Q_E^2} \right] \sum_{\{K_E\}} \frac{1}{K_E^2} = \frac{2}{(d-2)} \hat{\mathcal{I}}_1 \mathcal{I}_2, \quad (C.13) $$

$$ \hat{J}_5 = \hat{\mathcal{I}}_{1;111}^{01} = -\frac{5T^2}{1152\pi^2} \left[ \frac{1}{\epsilon} + 4\ln\frac{\bar{\mu}}{\pi T} + \frac{6}{5}\gamma_E + \frac{8}{3} + \frac{168}{5}\ln A + \frac{46}{5}\ln 2 - 0.199478 \right]. \quad (C.14) $$

We used the one-loop fermionic master integrals (C.9) and that $\hat{J}_2$ too vanishes identically, see e.g. [40]. Details on $\hat{J}_3$ and $\hat{J}_5$ are provided below in section C.1.1.

From eq. (3.22), eq. (B.4) and the definitions of the master integrals, we have

$$ Z_{g,q^{-2}}^{(c)} = -(D-2)g^2 \left[ N_c \left( (D-2)J_1 - J_2 - \frac{J_3}{2} - J_4 + 2J_5 \right) \right. \\ \left. \qquad - 2T_F N_F \left( 2\hat{J}_1 - \frac{D-4}{D-2}\hat{J}_2 - \frac{\hat{J}_3}{D-2} - \frac{2\hat{J}_4}{D-2} + \frac{4\hat{J}_5}{D-2} \right) \right]. \quad (C.15) $$

The bosonic part gives

$$ Z_{g,q^{-2},B}^{(c)} = \frac{5g^2 N_c T^2}{576\pi^2} \left[ -\frac{1}{\epsilon} - 4\ln\frac{\bar{\mu}}{4\pi T} + \frac{48}{5}\ln A - \frac{8}{5} \times [-0.299378] \right. \\ \left. \qquad - \frac{24}{5}\gamma_E - \frac{13}{15} - \frac{36}{5}\ln 2 \right], \quad (C.16) $$

$$ Z_{g,q^{-2},B,n=0}^{(c)} = \frac{3g^2 N_c T^2}{64\pi^2} \left[ -\frac{1}{\epsilon} - 4\ln\left(\frac{\bar{\mu}}{2T}\right) - \frac{5}{3} \right], \quad (C.17) $$

where the last line is the $Q_E$ zero-mode part of the $J_i$ integrals, before any $Q_E$ shift has been performed. It gives rise to the soft IR divergence and can be evaluated relatively straightforwardly (cf. e.g. eq. (C.44) and [41]).

After subtracting this zero-mode from the bosonic contribution — which we shall match with EQCD later — we find the bosonic and fermionic parts of eq. (C.15)

$$ \begin{aligned} Z_{g,q^{-2},B,n \neq 0}^{(c)} &= Z_{g,q^{-2},B}^{(c)} - Z_{g,q^{-2},B,n=0}^{(c)} \\ &= \frac{g^2 N_c T^2}{96\pi^2} \frac{11}{3} \left[ \frac{1}{\epsilon} + 4\ln\frac{\bar{\mu}}{T} + \frac{24}{11}\ln A - \frac{4}{11} \times [-0.299378] \right. \\ & \qquad \left. - \frac{12}{11}\gamma_E + \frac{61}{33} + \frac{10}{11}\ln\pi - \frac{52}{11}\ln 2 \right], \end{aligned} \quad (C.18) $$

$$ Z_{g,q^{-2},F}^{(c)} = -\frac{g^2 T_F N_F T^2}{72\pi^2} \left[ \frac{1}{\epsilon} + 4\ln\frac{\bar{\mu}A^3}{4\pi T} + \frac{5}{2} \times [-0.199478] + 3\gamma_E + \frac{13}{6} \right]. \quad (C.19) $$

### C.1.1 Evaluation of $J_3$ and $J_5$ and their fermionic counterparts

For the bosonic $J_3$ and fermionic $\hat{J}_3$, we first perform the Matsubara sums and drop the scale-free vacuum × vacuum term, i.e.

$$ J_3 = \sum_{Q_E K_E} \frac{1}{Q_E^2 K_E^2 (\mathbf{q} - \mathbf{k})^2} = \int_{\mathbf{q}\mathbf{k}} \frac{[1 + 2n_B(q)][1 + 2n_B(k)]}{4q\mathbf{k}(\mathbf{q} - \mathbf{k})^2} = \int_{\mathbf{q}\mathbf{k}} \frac{n_B(q)[1 + n_B(k)]}{q\mathbf{k}(\mathbf{q} - \mathbf{k})^2}, \quad (C.20) $$

$$ \hat{J}_3 = \sum_{\{Q_E K_E\}} \frac{1}{Q_E^2 K_E^2 (\mathbf{q} - \mathbf{k})^2} = - \int_{\mathbf{q}\mathbf{k}} \frac{n_F(q)[1 - n_F(k)]}{q\mathbf{k}(\mathbf{q} - \mathbf{k})^2}, \quad (C.21) $$
where we have used the $\mathbf{q} \leftrightarrow \mathbf{k}$ symmetry of the integrand. We can now split the vacuum and thermal parts in $k$. The former can be treated as follows

$$J_3^{(0)} = + \int_{\mathbf{q}\mathbf{k}} \left[ \frac{n_B(q)}{q\mathbf{k}(\mathbf{q}-\mathbf{k})^2} - \frac{n_B(q)\theta(k-q)}{q\mathbf{k}^3} + \frac{n_B(q)\theta(k-q)}{q\mathbf{k}^3} \right], \quad (C.22)$$

$$\hat{J}_3^{(0)} = - \int_{\mathbf{q}\mathbf{k}} \left[ \frac{n_F(q)}{q\mathbf{k}(\mathbf{q}-\mathbf{k})^2} - \frac{n_F(q)\theta(k-q)}{q\mathbf{k}^3} + \frac{n_F(q)\theta(k-q)}{q\mathbf{k}^3} \right]. \quad (C.23)$$

The sum of the first two terms in both cases is now finite and can be treated in $d=3$, whereas the third term is doable in DR. We then find

$$J_3^{(0)} = \frac{T^2}{24\pi^2} + \frac{T^2}{48\pi^2} \left[ \frac{1}{\epsilon} - 4 \ln \frac{4\pi T}{\bar{\mu} A^{12}} \right], \quad (C.24)$$

$$\hat{J}_3^{(0)} = -\frac{T^2}{48\pi^2} - \frac{T^2}{96\pi^2} \left[ \frac{1}{\epsilon} - 4 \ln \frac{8\pi T}{\bar{\mu} A^{12}} \right], \quad (C.25)$$

where the first term comes from the finite part and the second from the DR part.

We now turn to the thermal $k$ contribution. There we need to subtract its zero mode only for the bosonic integral and treat it separately in DR, i.e.

$$J_3^{(T)} = \int_{\mathbf{q}\mathbf{k}} \left[ \frac{n_B(q)n_B(k)}{q\mathbf{k}(\mathbf{q}-\mathbf{k})^2} - \frac{n_B(q)T}{q\mathbf{k}^2(\mathbf{q}-\mathbf{k})^2} + \frac{n_B(q)T}{q\mathbf{k}^2(\mathbf{q}-\mathbf{k})^2} \right], \quad \hat{J}_3^{(T)} = \int_{\mathbf{q}\mathbf{k}} \frac{n_F(q)n_F(k)}{q\mathbf{k}(\mathbf{q}-\mathbf{k})^2}. \quad (C.26)$$

Once again, the sum of the first two terms is finite. It can be treated numerically in three dimensions. The third term is dealt with in DR, using standard vacuum techniques for the $\mathbf{k}$ integration followed by the $\mathbf{q}$ one. We find

$$J_3^{(T)} = \frac{T^2}{64\pi^2} \left[ -3.25532 - \frac{1}{\epsilon} - 4 \ln \frac{\bar{\mu}}{T} - 2 \right]. \quad (C.27)$$

Using the position-space method of [40, 42], the numerical part can be integrated analytically with full agreement, finding

$$J_3^{(T)} = \frac{T^2}{48\pi^2} \left[ 1 + 3 \ln \frac{2\pi}{A^{12}} \right] + \frac{T^2}{64\pi^2} \left[ -\frac{1}{\epsilon} - 4 \ln \frac{\bar{\mu}}{T} - 2 \right], \quad (C.28)$$

$$\hat{J}_3^{(T)} = \frac{T^2}{48\pi^2} [1 - \ln 2]. \quad (C.29)$$

Hence

$$J_3 = J_3^{(0)} + J_3^{(T)} = \frac{T^2}{192\pi^2} \left[ \frac{1}{\epsilon} - 4 \ln \frac{32\pi T}{\bar{\mu} A^{12}} + 6 \right], \quad (C.30)$$

$$\hat{J}_3 = \hat{J}_3^{(0)} + \hat{J}_3^{(T)} = -\frac{T^2}{96\pi^2} \left[ \frac{1}{\epsilon} - 4 \ln \frac{8\pi T}{\bar{\mu} A^{12}} + \ln 4 \right]. \quad (C.31)$$

To evaluate $J_5$, we employ again the positon-space methods of [42] and write $J_5$ as

$$J_5 = \oint_{Q_E K_E} \frac{k_0^2}{Q_E^2 q^2 K_E^2 (Q_E + K_E)^2} = \oint_{Q_E} \frac{1}{Q_E^2 q^2} \Pi(Q_E), \quad (C.32)$$
where we can separate the vacuum and thermal contribution to $\Pi(Q)$

$$ \Pi(Q_E) \equiv \oint_{K_E} \frac{k_0^2}{K_E^2(Q_E + K_E)^2} = \Pi^{(0)}(Q_E) + \Pi^{(T)}(Q_E). \quad (C.33) $$

The vacuum contribution is obtainable by replacing the Matsubara sum over $k_0$ with a Euclidean $\int dk_0/(2\pi)$. Employing standard DR methods and inserting the result into eq. (C.32), yields

$$ J_5^{(0)} = \mu^{6-2d} \frac{((d-2)d-1)T^{2d-4}}{(4\pi)^{d+1}\pi^{\frac{7}{2}-2d}} \frac{\Gamma(\frac{5}{2}-d)\Gamma(d-3)}{\Gamma(\frac{d}{2})\Gamma(\frac{d+2}{2})} \zeta_{5-2d} = -\frac{T^2}{576\pi^2} \left[ \frac{1}{\epsilon} + 4\ln\frac{\mu A^{12}}{4\pi T} + \frac{2}{3} \right]. \quad (C.34) $$

For $\Pi^{(T)}$, we instead proceed in position space, i.e.

$$ \begin{aligned} \Pi(Q_E) &= \oint_{K_E} \frac{k_0^2}{K_E^2(Q_E + K_E)^2} = \int d^3\mathbf{r} \oint_{K_E P_E} \frac{k_0^2 \beta \delta_{p_0-k_0-q_0}}{K_E^2 P_E^2} e^{i\mathbf{r} \cdot (\mathbf{p}+\mathbf{q}+\mathbf{k})} \\ &= T \sum_{k_0} \int d^3\mathbf{r} \frac{k_0^2 e^{-|k_0|r} e^{-|q_0+k_0|r}}{(4\pi r)^2} e^{i\mathbf{q} \cdot \mathbf{r}}, \end{aligned} \quad (C.35) $$

where $\delta_{p_0} = \delta_{p_0,0}$ and in the last step we took the three-dimensional Fourier transforms. One can also work out the Matsubara sum, finding

$$ \Pi(Q_E) = T^3 \int d^3\mathbf{r} \frac{e^{-|\bar{q}_0|\bar{r}} e^{i\mathbf{q}\cdot\mathbf{r}}}{12r^2} \left[ |\bar{q}_0|^3 + \frac{3}{2} \left( \bar{q}_0^2 \coth(\bar{r}) + \text{csch}^2(\bar{r}) (|\bar{q}_0| + \coth(\bar{r})) \right) + \frac{|\bar{q}_0|}{2} \right], \quad (C.36) $$

where $\bar{q}_0 \equiv q_0/(2\pi T)$ and $\bar{r} \equiv 2\pi Tr$. To get $\Pi^{(T)}$, we must subtract the vacuum contribution, again obtained by replacing the sum with the integral. We then find

$$ \Pi^{(T)}(Q_E) = T^3 \int d^3\mathbf{r} \frac{e^{-|\bar{q}_0|\bar{r}} e^{i\mathbf{q}\cdot\mathbf{r}}}{24r^2} \left[ 3 \left( \bar{q}_0^2 \left( \coth(\bar{r}) - \frac{1}{\bar{r}} \right) + |\bar{q}_0| \left( \text{csch}^2(\bar{r}) - \frac{1}{\bar{r}^2} \right) \right) + \text{csch}^2(\bar{r}) \coth(\bar{r}) - \frac{1}{\bar{r}^3} \right] + |\bar{q}_0| \quad (C.37) $$

In principle, we should naively insert this into eq. (C.32). We would find that

1) Since the Fourier transform

$$ \int_{\mathbf{q}} \frac{e^{i\mathbf{q}\cdot\mathbf{r}}}{q^2 Q_E^2} = \frac{1 - e^{-|q_0|r}}{4\pi q_0^2 r} (1 - \delta_{q_0}), \quad (C.38) $$

is only defined in 3D for the non-zero modes, we need to treat the $Q_E$ zero mode separately.

2) If we use eqs. (C.37) and (C.38) together, perform the nonzero $q_0$ sum, we are then faced with a $dr/r$ UV logarithmic divergence. As per [42], we must then subtract off the $Q_E \gg T$ limit of eq. (C.37) in position space and add it back in momentum space in DR.
We then find

$$ \Pi^{(T)}(Q_E \gg T) = T^3 \int d^3\mathbf{r} \frac{e^{-|\bar{q}_0|\bar{r}} e^{i\mathbf{q}\cdot\mathbf{r}}}{24r^2} \bar{q}_0^2 \bar{r} + \dots, \qquad (C.39) $$

where the dots stand for higher-order terms in the expansion for $\bar{r} \ll 1$ with $\bar{r}\bar{q}_0 \sim 1$. In momentum space the analogue expansion can be made by having either of the $K_E$ and $Q_E + K_E$ denominators in eq. (C.32) be much larger than $T$. The leading term comes from the former, yielding

$$ \Pi^{(T)}(Q_E \gg T) = \frac{q_0^2}{Q_E^2} \mathcal{F}_{K_E} \left( \frac{1}{K_E^2} + \mathcal{O}\left(\frac{T^4}{Q_E^2}\right) \right). \qquad (C.40) $$

By splitting $J_5^{(T)}$ into three terms and noting that primed sum-integrals exclude the zero mode, we obtain

$$ \begin{aligned} J_5^{(T)} &= \mathcal{F}_{Q_E}' \left[ \frac{\Pi^{(T)}(Q_E)}{Q_E^2 q^2} - \frac{\Pi^{(T)}(Q_E \gg T)}{Q_E^2 q^2} \right] + \mathcal{F}_{Q_E}' \frac{q_0^2}{Q_E^4 q^2} \mathcal{F}_{K_E} \left( \frac{1}{K_E^2} + \mathcal{F}_{Q_E K_E} \frac{\delta_{q_0} k_0^2}{q^4 K_E^2 (Q_E + K_E)^2} \right) \\ &= J_5^{(T) \text{ finite}} + J_5^{(T) \text{ div.}} + J_5^{(T), q_0=0}. \end{aligned} \qquad (C.41) $$

For the first term we proceed in position space. After using eq. (C.38) and performing the $q_0$ sum, we have

$$ \begin{aligned} J_5^{(T) \text{ finite}} &= T^4 \sum_{q_0 \neq 0} \int d^3\mathbf{r} \frac{(1-e^{-|\bar{q}_0|\bar{r}})e^{-|\bar{q}_0|\bar{r}}}{4\pi q_0^2 r} \frac{1}{8r^2} \left[ \left( \text{csch}^2(\bar{r}) \coth(\bar{r}) - \frac{1}{\bar{r}^3} \right) \right. \\ & \qquad \left. + |\bar{q}_0| \left( \text{csch}^2(\bar{r}) - \frac{1}{\bar{r}^2} + \frac{1}{3} \right) + |\bar{q}_0|^2 \left( \coth(\bar{r}) - \frac{1}{\bar{r}} - \frac{\bar{r}}{3} \right) \right] \\ &= \frac{T^5}{8} \int d^3\mathbf{r} \left\{ \left[ \text{Li}_2(e^{-\bar{r}}) - \text{Li}_2(e^{-2\bar{r}}) \right] \left( \frac{\coth(\bar{r})\text{csch}^2(\bar{r})}{\bar{r}^3} - \frac{1}{\bar{r}^6} \right) \right. \\ & \qquad \left. + \ln(e^{-\bar{r}}+1) \left( \frac{\text{csch}^2(\bar{r})}{\bar{r}^3} + \frac{1}{3\bar{r}^3} - \frac{1}{\bar{r}^5} \right) + \frac{\text{csch}(\bar{r})}{6\bar{r}^2} \left( \frac{3\coth(\bar{r})}{\bar{r}} - 1 - \frac{3}{\bar{r}^2} \right) \right\} \\ &= [-0.0997928] \times \frac{T^2}{96\pi^2}, \end{aligned} \qquad (C.42) $$

where the integration is finite and has been carried out numerically.<fnref n="7" />

For the divergent piece we find, using IBP

$$ \begin{aligned} J_5^{(T) \text{ div.}} &= \mathcal{F}_{Q_E}' \frac{q_0^2}{Q_E^4 q^2} \mathcal{F}_{K_E} \left( \frac{1}{K_E^2} \right) = \left[ \frac{4-d}{2} \mathcal{F}_{Q_E} \left( \frac{1}{Q_E^2 q^2} \right) \right] \mathcal{I}_1 = -\frac{d-4}{d-2} \mathcal{I}_1 \mathcal{I}_2 \\ &= \frac{T^2}{192\pi^2} \left[ \frac{1}{\epsilon} + 4 \ln \frac{\bar{\mu}A^6}{4\pi T} + 2\gamma_E + 4 \right]. \end{aligned} \qquad (C.43) $$

$$ \frac{T^5}{8} \int d^3\mathbf{r} \frac{\text{csch}(\bar{r})}{6\bar{r}^2} \left( \frac{3 \coth(\bar{r})}{\bar{r}} - 1 - \frac{3}{\bar{r}^2} \right) = \frac{T^2}{96\pi^2} \left[ \ln \frac{A^{12}}{8} - \frac{1}{2} - \gamma_E \right]. $$

<footnote n="7">⁷An analytical evaluation could be possible by using IBP: $\ln(1+e^{-\bar{r}}) = \ln(1-e^{-2\bar{r}}) - \ln(1-e^{-\bar{r}})$, and these two terms are related to the derivatives of the polylogarithms. The non-logarithmic and non-polylogarithmic term can be done using the methods of [40, 42], giving</footnote>
Finally, the zero mode can be computed using Feynman parametrisation and the series expansion for the Riemann zeta function. Below, we list the bosonic [43] and fermionic results

$$ \begin{aligned} \mathcal{Z}_{s_1 s_2 s_3}^{\alpha} &= \oint_{Q_E K_E} \delta_{q_0} \frac{(k_0^2)^{\alpha}}{[Q_E^2]^{s_1} [K_E^2]^{s_2} [(Q_E - K_E)^2]^{s_3}} \\ &= \mu^{6-2d} \frac{\Gamma(\frac{d}{2} - s_1) \Gamma(s_12 - \frac{d}{2}) \Gamma(s_13 - \frac{d}{2}) \Gamma(s_123 - d)}{(4\pi)^d \Gamma(\frac{d}{2}) \Gamma(s_2) \Gamma(s_3) \Gamma(s_{1123} - d)} \frac{2T^2 \zeta_{2s_{123}-2\alpha-2d}}{(2\pi T)^{2s_{123}-2\alpha-2d}}, \end{aligned} \quad (C.44) $$

$$ \hat{\mathcal{Z}}_{s_1 s_2 s_3}^{\alpha} = \oint_{Q_E \{K_E\}} \delta_{q_0} \frac{(k_0^2)^{\alpha}}{[Q_E^2]^{s_1} [K_E^2]^{s_2} [(Q_E - K_E)^2]^{s_3}} = (2^{2s_{123}-2\alpha-2d} - 1) \mathcal{Z}_{s_1 s_2 s_3}^{\alpha}, \quad (C.45) $$

where $s_{\{i\}} = \sum_{j \in \{i\}} s_j$. The two special cases relevant for our computation

$$ \mathcal{Z}_{211}^{1} \stackrel{d=3-2\epsilon}{\equiv} J_{5}^{(T)} q_{0}=0=\frac{T^{2}}{128 \pi^{2}}, \quad (C.46) $$

$$ \hat{\mathcal{Z}}_{211}^{1} \stackrel{d=3-2\epsilon}{\equiv} \hat{J}_{5}^{(T)} q_{0}=0=0, \quad (C.47) $$

can also be read from eqs. (D2) and (D5) of [41]. Summing up eqs. (C.34), (C.42), (C.43) and (C.46) we find eq. (C.7).

For $\hat{J}_5$ we adopt the same strategy used for its bosonic counterpart. The $K_E$ vacuum part is identical and given by eq. (C.34). The $K_E$ thermal part is split in three terms (cf. eq. (C.41))

$$ \hat{J}_{5}^{(T)} = \hat{J}_{5}^{(T) \text{ finite}} + \hat{J}_{5}^{(T) \text{ div.}} + \hat{J}_{5}^{(T), q_{0}=0}. \quad (C.48) $$

The finite and divergent terms are

$$ \begin{aligned} \hat{J}_5^{(T) \text{ finite}} &= \frac{T^5}{8} \int d^3 \mathbf{r} \left\{ \left[ \text{Li}_2(e^{-\bar{r}}) - \text{Li}_2(e^{-2\bar{r}}) \right] \left( \frac{\text{csch}^3(\bar{r})(\cosh(2\bar{r}) + 3)}{4\bar{r}^3} - \frac{1}{\bar{r}^6} \right) \right. \\ & \qquad \left. + \ln(e^{-\bar{r}} + 1) \left( \frac{6 \coth(\bar{r}) \cosh(\bar{r}) - 1}{6\bar{r}^3} - \frac{1}{\bar{r}^5} \right) \right. \\ & \qquad \left. + \frac{\cosh(\bar{r})}{12\bar{r}^2} \left( \frac{6 \cosh(\bar{r})}{\bar{r}} + 1 - \frac{6}{\bar{r}^2} \right) \right\} \\ &= [0.0831159] \times \frac{T^2}{96\pi^2}, \end{aligned} \quad (C.49) $$

$$ \hat{J}_5^{(T) \text{ div.}} = (2^{2-d} - 1) J_5^{(T) \text{ div.}} \stackrel{(C.43)}{=} - \frac{T^2}{384 \pi^2} \left[ \frac{1}{\epsilon} + 4 \ln \frac{\mu A^6}{4\pi T} + 2\gamma_E + 4 - \ln 4 \right]. \quad (C.50) $$

Finally, the zero mode can be read off from eq. (C.47). Summing up eqs. (C.34), (C.49), (C.50) and (C.47) we find eq. (C.14).

## C.2 Real-time evaluation of $Z_{g,L-T}$

For the sake of readability, we split our evaluation into the fermionic quark- and the bosonic gauge-contribution.
### C.2.1 Quark contribution

The unresummed Hard Thermal Loop (HTL) subset of the full contribution, using eq. (B.6), is

$$Z_{g,L-T,HTL,F}^{(c)} = i \int_Q \frac{1+n_B(q^0)}{(q^--i\epsilon)^2} \left[ \frac{q_{\perp}^2 (\Pi_L^R(Q)-\Pi_T^R(Q))_{HTL}}{q^2(Q^2+i\epsilon q^0)} - \text{adv.} \right], \quad (C.51)$$

$$\Pi_L^R(Q)_{\text{F,HTL}} - \Pi_T^R(Q)_{\text{F,HTL}} = \frac{g^2 T_F N_F T^2}{12} \left[ -\frac{(6q_0^2 - 4q^2)}{q^2} + 3 \frac{q^0 Q^2}{q^3} \ln \frac{q^0 + q + i\epsilon}{q^0 - q + i\epsilon} \right], \quad (C.52)$$

which, after inserting eq. (C.52) into (C.51), gives rise to a pole and a cut contribution. The pole contribution is

$$Z_{g,L-T,F,HTL \text{ pole}}^{(c)} = -\frac{g^2 T_F N_F T^2}{12} \int_Q \frac{\frac{1}{2} + n_B(|q^0|)}{(q^--i\epsilon)^2} \frac{(6q_0^2 - 4q^2)}{q^2} \frac{q_{\perp}^2}{q^2} 2\pi\delta(Q^2) \quad (C.53)$$

$$= -\frac{g^2 T_F N_F T^2}{6} \int_q \frac{n_B(q)}{2q} \left[ \frac{1}{(q-q^z-i\epsilon)^2} + \frac{1}{(-q-q^z-i\epsilon)^2} \right] \frac{q_{\perp}^2}{q^2} \quad (C.54)$$

$$= \frac{g^2 T_F N_F T^2}{48\pi^2} \left[ \frac{1}{\epsilon^2} + \frac{2}{\epsilon} \ln \frac{\bar{\mu} e^{\gamma_E + \frac{1}{2}}}{4\pi T} + 2 \ln^2 \frac{\bar{\mu} e^{\gamma_E + \frac{1}{2}}}{4\pi T} - 4\gamma_1 + \frac{\pi^2}{4} - 2\gamma_E^2 + \frac{3}{2} \right],$$

where in DR, the second step used that the vacuum part is scale-free. The cut contribution of the HTL piece is instead

$$Z_{g,L-T,F,HTL \text{ cut}}^{(c)} = \frac{g^2 T_F N_F T^2}{4} \int_q \int_0^q \frac{dq^0}{2\pi} \frac{q_{\perp}^2 q^0}{q^5} 2\pi \left[ \frac{1+n_B(q^0)}{(q^0-q^z-i\epsilon)^2} + \frac{n_B(q^0)}{(q^0+q^z+i\epsilon)^2} \right], \quad (C.55)$$

where we have symmetrised the frequency integral before restricting the integration to be over twice the positive range. The $n_B$-independent part is now scale-free, while the rest can be done in a somewhat similar way to the pole piece

$$Z_{g,L-T,F,HTL \text{ cut}}^{(c)} = \frac{g^2 T_F N_F T^2}{4} \left(\frac{\bar{\mu}^2 e^{\gamma_E}}{4\pi}\right)^\epsilon \frac{2\pi^{\frac{d-1}{2}}}{\Gamma(\frac{d-1}{2})(2\pi)^{d-1}} \int_0^\infty \frac{dq^0}{2\pi} q^0 n_B(q^0) \\ \times \int_{q_0}^\infty dq q^{d-4} \int_{-1}^1 dx (1-x^2)^{\frac{d-1}{2}} \left[ \frac{1}{(q^0-qx-i\epsilon)^2} + \frac{1}{(q^0+qx+i\epsilon)^2} \right]. \quad (C.56)$$

Our strategy here is to do the following subtraction: $n_B(q^0) \to n_B(q^0) \mp (T/q^0 - \theta(T - q^0)/2)$. In this way, the subtraction is finite, and then we re-evaluate the addition in DR. For the finite part, we find, performing the angular and momentum integrations first

$$Z_{g,L-T,F,HTL \text{ cut finite}}^{(c)} = \frac{g^2 T_F N_F T^2}{4} \frac{1}{2\pi} \int_0^\infty \frac{dq^0}{2\pi} q^0 \left[ n_B(q^0) - \frac{T}{q^0} + \frac{\theta(T-q^0)}{2} \right] \\ \times \int_{q^0}^\infty \frac{dq}{q} \int_{-1}^1 dx (1-x^2) \left[ \frac{1}{(q^0-qx-i\epsilon)^2} + \frac{1}{(q^0+qx+i\epsilon)^2} \right] \\ = \frac{g^2 T_F N_F T^2}{12\pi^2} \ln \frac{e^{\gamma_E}}{2\pi} (\ln 2 - 1). \quad (C.57)$$
For the DR subtracted part, we instead do the frequency integral first, splitting it into the $q > T$ and $q < T$ ranges,

$$Z_{g,L-T,F,HTL \text{ cut div}}^{(c)} = -\frac{g^2 T_F N_F T^2}{8} \left(\frac{\bar{\mu}^2 e^{\gamma_E}}{4\pi}\right)^\epsilon \frac{2\pi^{\frac{d-1}{2}}}{\Gamma\left(\frac{d-1}{2}\right) (2\pi)^d} \int_0^\infty dq \, q^{d-4} \int_{-1}^1 dx (1-x^2)^{\frac{d-1}{2}} \\ \times \left\{ \theta(T-q) \left[ -\mathbb{P}\frac{2}{1-x^2} + \ln\frac{1-x^2}{x^2} \right] \right. \\ \left. + \theta(q-T) \left[ \ln\frac{|T^2-q^2x^2|}{q^2x^2} - \mathbb{P}\frac{2T^2}{T^2-q^2x^2} \right] \right\}, \quad (C.58)$$

where the principal value (PV) $\mathbb{P}$ arises from the $\epsilon \to 0$ limit. The $q > T$ slice is finite whereas the $T > q$ part is IR-divergent and needs regularisation

$$Z_{g,L-T,F,HTL \text{ cut div}}^{(c)} = -\frac{g^2 T_F N_F T^2}{8} \left(\frac{\bar{\mu}^2 e^{\gamma_E}}{4\pi}\right)^\epsilon \frac{2\pi^{\frac{d-1}{2}}}{\Gamma\left(\frac{d-1}{2}\right) (2\pi)^d} \\ \times \left\{ \int_{-1}^{1} dx (1-x^2)^{\frac{d-1}{2}} \frac{T^{d-3}}{d-3} \left[ -\mathbb{P}\frac{2}{1-x^2} + \ln\frac{1-x^2}{x^2} \right] \right. \\ \left. + \int_{-1}^{1} dx (1-x^2) \left[ \frac{\text{Li}_2(x^2)}{2} + \ln\left(\frac{x^2}{1-x^2}\right) + \ln^2|x| - \frac{\pi^2}{6} \right] \right\}, \quad (C.59) \\ = -\frac{g^2 T_F N_F T^2}{8} \left(\frac{\bar{\mu}^2 e^{\gamma_E}}{4\pi}\right)^\epsilon \frac{2\pi^{\frac{d-1}{2}}}{\Gamma\left(\frac{d-1}{2}\right) (2\pi)^d} \\ \times \left\{ \frac{T^{d-3}}{d-3} \frac{\Gamma(\frac{1}{2})\Gamma(\frac{d+1}{2})}{\Gamma(\frac{d}{2}+1)} \left( \psi(\frac{d-1}{2}) + \gamma_E - 2 + 2\ln 2 \right) - \frac{1}{9} (8+\pi^2-8\ln 2) \right\} \\ = \frac{g^2 T_F N_F T^2}{24\pi^2} (\ln 2 - 1) \left[ \frac{1}{\epsilon} + 2\ln\frac{\bar{\mu}}{2T} + 1 \right], \quad (C.60)$$

where $\psi(x) = (\ln \Gamma(x))'$ is the digamma function. Finally

$$Z_{g,L-T,F,HTL \text{ cut}}^{(c)} = \frac{g^2 T_F N_F T^2}{24\pi^2} (\ln 2 - 1) \left[ \frac{1}{\epsilon} + 2 \ln \frac{\bar{\mu} e^{\gamma_E + \frac{1}{2}}}{4\pi T} \right]. \quad (C.61)$$

Let us now consider the non-HTL piece

$$\delta\Pi_L^R(Q)_F - \delta\Pi_T^R(Q)_F = \frac{2g^2 T_F N_F Q^2}{(4\pi)^2 q^3} \int_0^\infty dk \, n_F(k) \left[ 24kq \right. \\ \left. + (3(2k-q^0)^2 - q^2) \ln \frac{q^0 + q - 2k + i\epsilon}{q^0 - q - 2k + i\epsilon} \right. \\ \left. + (3(2k+q^0)^2 - q^2) \ln \frac{q^0 - q + 2k + i\epsilon}{q^0 + q + 2k + i\epsilon} \right], \quad (C.62)$$
where $\delta\Pi = \Pi - \Pi_{HTL}$. This then gives

$$Z_{g,L-T,F,\text{non-HTL}}^{(c)} = \frac{g^2 T_F N_F}{4\pi} \int_0^\infty dk \, n_F(k) \int_Q \frac{1}{(q^- - i\epsilon)^2} \frac{q_\perp^2}{q^5} \\ \times \left[ \theta(q^2 - (q^0 - 2k)^2) (3(2k - q^0)^2 - q^2) \right. \\ \left. - \theta(q^2 - (q^0 + 2k)^2) (3(2k + q^0)^2 - q^2) \right]. \qquad (C.63)$$

Here, we used the fact that the $\frac{1}{2}$ piece vanishes$^8$ to restrict to the purely odd $\frac{1}{2} + n_B(q^0)$ component. Let us again symmetrise and restrict to twice the positive frequency range

$$Z_{g,L-T,F,\text{non-HTL}}^{(c)} = \frac{g^2 T_F N_F}{4\pi} \int_0^\infty dk \, n_F(k) \int_0^\infty \frac{dq^0}{2\pi} \int_{\mathbf{q}} \frac{q_\perp^2}{q^5} \left[ \frac{1}{(q^0 - q^z - i\epsilon)^2} + \frac{1}{(q^0 + q^z + i\epsilon)^2} \right] \\ \times \left\{ \theta(q-2k) \left[ (3(2k-q^0)^2-q^2)\theta(2k+q-q^0) \right. \right. \\ \left. \left. - (3(2k+q^0)^2-q^2)\theta(q-2k-q^0) \right] \right. \\ \left. + \theta(2k-q) (3(2k-q^0)^2-q^2)\theta(2k+q-q^0)\theta(q^0+q-2k) \right\}. \qquad (C.64)$$

We start by dealing with the vacuum part

$$Z_{g,L-T,F,\text{non-HTL vac}}^{(c)} = \frac{g^2 T_F N_F}{8\pi} \int_0^\infty dk \, n_F(k) \int_0^\infty \frac{dq^0}{2\pi} \int_{\mathbf{q}} \frac{q_\perp^2}{q^5} \left[ \frac{1}{(q^0 - q^z - i\epsilon)^2} + \frac{1}{(q^0 + q^z + i\epsilon)^2} \right] \\ \times \left\{ \theta(q-2k) \left[ (3(2k-q^0)^2-q^2)\theta(q+2k-q^0) \right. \right. \\ \left. \left. - (3(2k+q^0)^2-q^2)\theta(q-2k-q^0) \right] \right. \\ \left. + \theta(2k-q) (3(2k-q^0)^2-q^2)\theta(2k+q-q^0)\theta(q^0-2k+q) \right\}. \qquad (C.65)$$

The strategy here is the following: we further symmetrise the expression in $q^z$, to perform the $q^0$ integrations. After that one can take $\epsilon \to 0$. Logarithmic pieces will see their argument change to its absolute value, while possible imaginary parts will cancel out due to $q^z$ symmetrisation. Non-logarithmic pieces present poles on the $q^z$ integration range. If we restrict to twice the positive $q^z$ range, these occur at $-|2k - q| + q^z = 0$ and are turned into a PV prescription by the $\epsilon \to 0$ limit. We can thus proceed to the $q^z$ integral in the two ranges $2k > q$ and $q > 2k$. In the latter case, we find two subranges, one for $q_\perp < 2k$ and one for $q_\perp > 2k$. This last one has to be carried out in DR, since it contains a UV divergence. We then find

$$Z_{g,L-T,F,\text{non-HTL vac}}^{(c)} = \frac{g^2 T_F N_F T^2}{48\pi^2} \left\{ \frac{1}{\epsilon^2} + \frac{2}{\epsilon} \ln \frac{\bar{\mu} A^{12} e^{-\frac{1}{2}}}{16\pi T} + 2 \ln^2 \frac{\bar{\mu} A^{12} e^{-\frac{1}{2}}}{16\pi T} - \frac{7}{2} \ln^2(2) - 2(\ln \zeta_2)^2 \right. \\ \left. + \frac{2\zeta_2''}{\zeta_2} + \frac{1}{2} \text{Li}_2\left(\frac{2}{3}\right) + \frac{1}{4} \text{Li}_2\left(\frac{3}{4}\right) - \frac{1}{2} - \frac{1}{4} \ln(3) \ln \frac{4}{3} \right\}. \qquad (C.66)$$

<footnote>$^8$It gives rise to fully retarded and fully advanced contributions in $q^+$. The contours can be safely closed in these half-planes without encountering singularities or contributions from the arcs at large $|q^+|$, thus yielding a vanishing contribution.</footnote>
For the thermal part, we have instead

$$
Z_{g,L-T,F,\text{non-HTL th}}^{(c)} = \frac{g^2 T_F N_F}{2\pi} \int_0^\infty dk \, n_F(k) \int_0^\infty \frac{dq^0}{2\pi} \int_0^\infty \frac{dq}{4\pi^2} \frac{n_B(q^0)}{q} \left[ \frac{2q^0 \tanh^{-1}\left(\frac{2qq^0}{q^2+q_0^2}\right)}{q^3} - \frac{4}{q^2} \right] \\
\times \left\{ \theta(q-2k) \left[ (3(2k-q^0)^2-q^2)\theta(2k+q-q^0) \right. \right. \\
\left. \left. - (3(2k+q^0)^2-q^2)\theta(q-2k-q^0) \right] \right. \\
\left. + \theta(2k-q) (3(2k-q^0)^2-q^2)\theta(2k+q-q^0)\theta(q^0+q-2k) \right\}. \quad (C.67)
$$

Here, we first perform the *q* integration analytically, followed by *q*⁰ and *k* numerically. We find

$$
Z_{g,L-T,F,\text{non-HTL th}}^{(c)} = \frac{g^2 T_F N_F T^2}{48\pi^2} \left[ \underbrace{-0.438129}_{q<2k} + \underbrace{1.97382}_{q>2k} \right]. \quad (C.68)
$$

Putting everything together, we find

$$
Z_{g,L-T,F}^{(c)} = \frac{g^2 T_F N_F T^2}{48\pi^2} \left\{ \begin{aligned}
& \frac{2}{\epsilon^2} + \frac{2}{\epsilon} \left( \ln \frac{\bar{\mu} e^{\gamma_E + \frac{1}{2}}}{4\pi T} + \ln \frac{\bar{\mu} A^{12} e^{-\frac{1}{2}}}{16\pi T} + \ln 2 - 1 \right) \\
& + 2 \ln^2 \frac{\bar{\mu} e^{\gamma_E + \frac{1}{2}}}{4\pi T} - 4\gamma_1 + \frac{\pi^2}{4} - 2\gamma_E^2 + \frac{3}{2} + 2(\ln 2 - 1) \left[ 2 \ln \frac{\bar{\mu} e^{\gamma_E + \frac{1}{2}}}{4\pi T} \right] \\
& + 2 \ln^2 \frac{\bar{\mu} A^{12} e^{-\frac{1}{2}}}{16\pi T} - \frac{7}{2} \ln^2(2) - 2(\ln \zeta_2)'^2 + \frac{2\zeta_2''}{\zeta_2} + \frac{1}{2} \text{Li}_2\left(\frac{2}{3}\right) + \frac{1}{4} \text{Li}_2\left(\frac{3}{4}\right) - \frac{1}{2} \\
& - \frac{1}{4} \ln(3) \ln \frac{4}{3} \underbrace{-0.438129}_{q<2k} + \underbrace{1.97382}_{q>2k}
\end{aligned} \right\}. \quad (C.69)
$$

C.2.2 Gluon and ghost contribution

We start from eqs. (3.20) and (B.5)

$$
Z_{g,L-T,B}^{(c)} = i \frac{g^2 C_A}{(4\pi)^2} \int_Q \int_0^\infty dk \frac{(1+n_B(q^0))}{q^3(q^--i\epsilon)^2} n_B(k) \left[ \frac{q_\perp^2}{q^2(Q^2+i\epsilon q^0)} \left( -8kq(3q_0^2-q^2) \right. \right. \\
\left. \left. +Q^2(12k^2-q^2+3q_0^2) \left( \ln \frac{q_0+q-2k+i\epsilon}{q_0-q-2k+i\epsilon} + \ln \frac{q_0-q+2k+i\epsilon}{q_0+q+2k+i\epsilon} \right) \right. \right. \\
\left. \left. +12kq_0Q^2 \left( 2\ln \frac{q_0+q+i\epsilon}{q_0-q+i\epsilon} + \ln \frac{q_0-q+2k+i\epsilon}{q_0+q+2k+i\epsilon} - \ln \frac{q_0+q-2k+i\epsilon}{q_0-q-2k+i\epsilon} \right) \right) - \text{adv.} \right]. \quad (C.70)
$$

As was done in section C.1, we need to separate out the contributions from zero and non-zero
Q Matsubara modes. In this case, because of the presence of the 1/(q⁻ - iε)² factor, the q⁰
integrand possesses poles in addition to those associated with the Matsubara frequencies,
meaning that we cannot effectively trade the q⁰ integral for a sum over Matsubara modes as
was done before. However, by first shifting: q² → q² + x/z²q⁰, we may eliminate the q⁰ pole
tied to the 1/(q⁻ - iε)² prefactor, before proceeding to carry out the q⁰ integration using
contour methods. Of course, one may worry that such a shift spoils the analytical structure
of the retarded (advanced) functions in the upper (lower) $q^0$ plane. It has nevertheless been pointed out by Caron-Huot [15]$^9$ that as long as $X$ is lightlike or spacelike, one may always work in a frame where $|\frac{x^0}{x^z}| \le 1$, in which case the analytical structure of the causal functions is not altered.$^{10}$ Thus, in this case where $X$ is lightlike, the appropriate shift is precisely: $q^z \to q^z + q^0$. Following this maneuver, we can safely carry out the integral $q^0$ and extract the $n=0$ contribution

$$ 
\begin{aligned}
Z_{g,L-T,B,n=0}^{(c)} &= \frac{g^2 C_A T}{8\pi^2} \int_{\mathbf{q}} \int_0^{\infty} dk \frac{q_{\perp}^2 n_B(k)}{q^5(q^z + i\epsilon)^2} \left[ 4kq - (12k^2 - q^2) \ln \left| \frac{2k-q}{q+2k} \right| \right] \\
&= \frac{g^2 C_A T^2}{32\pi^2} \left[ -\frac{1}{\epsilon} + 4 - \ln \frac{\bar{\mu}^2 e^{\gamma_E}}{4T^2} - \psi\left(\frac{3}{2}\right) \right],
\end{aligned}
\quad (C.71) $$

which is needed for the matching calculation completed in section 4.

For the contribution from the non-zero Matsubara modes, we rather persist in the Minkowskian signature so that we can employ the same strategy as in appendix C.2.1. In more detail, for the HTL part, the $k$ integrand is given by $kn_B(k)$. Therefore, given that

$$ \int_0^\infty dk \, k \, n_B(k) = 2 \int_0^\infty dk \, k \, n_F(k) = \frac{\pi^2 T^2}{6}, \quad (C.72) $$

we may simply make the replacement $2N_F T_F \to C_A$ in eq. (C.61) to obtain the bosonic equivalent. This yields

$$ Z_{g,L-T,B,HTL}^{(c)} = \frac{g^2 C_A T^2}{48\pi^2} \left[ \frac{1}{\epsilon^2} + \frac{2}{\epsilon} \ln \frac{\bar{\mu} e^{\gamma_E + \frac{1}{2}}}{4\pi T} + 2 \ln^2 \frac{\bar{\mu} e^{\gamma_E + \frac{1}{2}}}{4\pi T} - 4\gamma_1 + \frac{\pi^2}{4} - 2\gamma_E^2 + \frac{3}{2} \right] \\ + \frac{g^2 C_A T^2}{24\pi^2} (\ln 2 - 1) \left[ \frac{1}{\epsilon} + 2 \ln \frac{\bar{\mu} e^{\gamma_E + \frac{1}{2}}}{4\pi T} \right]. \quad (C.73) $$

For the vacuum part of the non-HTL contribution, we have instead

$$ Z_{g,L-T,B,\text{non-HTL vac}}^{(c)} = \frac{g^2 C_A}{16\pi} \int_0^\infty dk \, n_B(k) \int_0^\infty \frac{dq^0}{2\pi} \int_{\mathbf{q}} \frac{q_\perp^2}{q^5} \left[ \frac{1}{(q^0 - q^z - i\epsilon)^2} + \frac{1}{(q^0 + q^z + i\epsilon)^2} \right] \\ \times \left\{ \theta(q-2k) \left[ (3(2k-q^0)^2 - q^2) \theta(q+2k-q^0) \theta(q^0-q+2k) \right. \right. \\ \left. \left. - 24kq^0\theta(q-2k-q^0) \right] \right. \\ \left. + \theta(2k-q) (3(2k-q^0)^2 - q^2) \theta(2k+q-q^0) \theta(q^0-2k+q) \right\}, \quad (C.74) $$

$$ Z_{g,L-T,B,\text{non-HTL vac}}^{(c)} = \frac{g^2 C_A T^2}{48\pi^2} \left\{ \frac{1}{\epsilon^2} + \frac{2}{\epsilon} \ln \frac{\bar{\mu} A^{12} e^{-\frac{1}{2}}}{8\pi T} + 2 \ln^2 \frac{\bar{\mu} A^{12} e^{-\frac{1}{2}}}{8\pi T} + \frac{1}{2} \ln^2(2) - 2(\ln \zeta_2)^2 \right. \\ \left. + \frac{2\zeta_2''}{\zeta_2} + \frac{1}{2} \text{Li}_2\left(\frac{2}{3}\right) + \frac{1}{4} \text{Li}_2\left(\frac{3}{4}\right) - \frac{1}{2} - \frac{1}{4} \ln(3) \ln \frac{4}{3} \right\}, \quad (C.75) $$

<footnote>$^9$See also [44, 45].</footnote>

<footnote>$^{10}$While the methods introduced in [15] are strictly speaking only defined for $|\frac{x^0}{x^z}| < 1$, they can be safely extended to null correlators as well, provided the zero mode contribution is free of collinear singularities, which is the case here.</footnote>
where in the last step, the methods of appendix C.2.1 were used. For the thermal part, we have to subtract the zero mode (cf. eq. (C.71)$^{11}$). This yields

$$Z_{g,L-T,B,\text{non-HTL th},n\neq0}^{(c)} = \frac{g^2 C_A}{4\pi} \int_0^\infty dk \, n_B(k) \int_0^\infty \frac{dq^0}{2\pi} \int_0^\infty \frac{dq}{4\pi^2 q} \left[n_B(q^0) - \frac{T}{q^0}\right] \\ \times \left[\frac{2q^0 \tanh^{-1}\left(\frac{2qq^0}{q^2+q_0^2}\right)}{q^3} - \frac{4}{q^2}\right] \\ \times \left\{\theta(q-2k) \left[(3(2k-q^0)^2 - q^2)\theta(2k+q-q^0) - (3(2k+q^0)^2 - q^2)\theta(q-2k-q^0)\right] + \theta(2k-q) (3(2k-q^0)^2 - q^2)\theta(2k+q-q^0)\theta(q^0+q-2k)\right\}. \quad (C.76)$$

Here, we first perform the $q$ integral analytically, followed by the $q^0$ and $k$ integrals numerically. We find

$$Z_{g,L-T,B,\text{non-HTL th} n \neq 0}^{(c)} = \frac{g^2 C_A T^2}{48 \pi^2} \left[ \underbrace{0.145277}_{q<2k} - \underbrace{-5.13832}_{q>2k} \right]. \quad (C.77)$$

Not including the $Q$ zero-mode contribution, the total is then

$$Z_{g,L-T,B,n \neq 0}^{(c)} = \frac{g^2 C_A T^2}{48 \pi^2} \left\{ \frac{2}{\epsilon^2} + \frac{2}{\epsilon} \left( \ln \frac{\bar{\mu} e^{\gamma_E + \frac{1}{2}}}{4\pi T} + \ln \frac{\bar{\mu} A^{12} e^{-\frac{1}{2}}}{8\pi T} + \ln 2 - 1 \right) \right. \\ + 2 \ln^2 \frac{\bar{\mu} e^{\gamma_E + \frac{1}{2}}}{4\pi T} - 4\gamma_1 + \frac{\pi^2}{4} - 2\gamma_E^2 + \frac{3}{2} + 2(\ln 2 - 1) \left[ 2 \ln \frac{\bar{\mu} e^{\gamma_E + \frac{1}{2}}}{4\pi T} \right] \\ + 2 \ln^2 \frac{\bar{\mu} A^{12} e^{-\frac{1}{2}}}{8\pi T} + \frac{1}{2} \ln^2(2) - 2(\ln \zeta_2)'^2 + \frac{2\zeta_2''}{\zeta_2} + \frac{1}{2} \mathrm{Li}_2\left(\frac{2}{3}\right) + \frac{1}{4} \mathrm{Li}_2\left(\frac{3}{4}\right) - \frac{1}{2} \\ \left. - \frac{1}{4} \ln(3) \ln \frac{4}{3} + \underbrace{0.145277}_{q<2k} - \underbrace{-5.13832}_{q>2k} \right\}. \quad (C.78)$$

## D Determination of the QCD coupling

When adding perturbative and lattice determinations we need to be consistent in our choice of the coupling $g$. Indeed, we need to set it to the values that correspond to those chosen in the EQCD matching coefficients in the lattice ensembles. The matching procedure and coupling determination for these ensembles is detailed in [21, 22], following [20, 46]. The resulting values of the coupling and the number of light flavors $N_F$ are listed in table 2, together with the data for $Z_g^{3d}|_{\text{merge}}$ from [6].$^{12}$ The matching procedure corresponds to a renormalisation scale

$$\bar{\mu}_{UV} = 4\pi T e^{-\gamma_E} \hat{\mu}, \quad (D.1)$$

<footnote>$^{11}$No zero-mode subtraction was necessary in the HTL part, since the zero mode is scale-free there.</footnote>

<footnote>$^{12}$In the notation and normalisation of [6], that quantity is $(Z_g^{3d}/g_{3d})_{\text{full}}$. It is listed in table 1 there; please refer to the latest arXiv version, which fixes a numerical issue in previous versions.</footnote>
<table><caption>**Table 2.** The values for the QCD coupling squared and the number of light flavors correspond to those used in the matching coefficients in the lattice ensembles of [5, 6], see [21, 22] for details. Data for $Z_g^{3d}$'merge comes from [6].</caption><thead><tr><th>T</th><th>g<sup>2</sup></th><th>N<sub>F</sub></th><th>(Z<sub>g</sub><sup>3d</sup>)<sup>merge</sup> / (g<sub>3d</sub><sup>2</sup>T)</th></tr></thead><tbody><tr><td>250 MeV</td><td>3.725027</td><td>3</td><td>0.085(37)(12)(2)</td></tr><tr><td>500 MeV</td><td>2.763516</td><td>3</td><td>-0.024(36)(14)(1)</td></tr><tr><td>1 GeV</td><td>2.210169</td><td>4</td><td>-0.026(32)(4)(3)</td></tr><tr><td>100 GeV</td><td>1.066561</td><td>5</td><td>-0.179(15)(5)(2)</td></tr></tbody></table>

where at one-loop level [20]

$$ \hat{\mu} = \exp \frac{-N_c + 4 \ln(4) N_F}{22 N_c - 4 N_F}. \qquad (D.2) $$

Hence, we set $\bar{\mu}_{UV}$ to eq. (D.1) when numerically evaluating eq. (3.28) in section 5.

**Open Access.** This article is distributed under the terms of the Creative Commons Attribution License (CC-BY4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.

## References

[1] M. Connors, C. Nattrass, R. Reed and S. Salur, Jet measurements in heavy ion physics, Rev. Mod. Phys. 90 (2018) 025005 [arXiv:1705.01974] [INSPIRE].

[2] L. Cunqueiro and A.M. Sickles, Studying the QGP with Jets at the LHC and RHIC, Prog. Part. Nucl. Phys. 124 (2022) 103940 [arXiv:2110.14490] [INSPIRE].

[3] L. Apolinário, Y.-J. Lee and M. Winn, Heavy quarks and jets as probes of the QGP, Prog. Part. Nucl. Phys. 127 (2022) 103990 [arXiv:2203.16352] [INSPIRE].

[4] P.B. Arnold, G.D. Moore and L.G. Yaffe, Photon and gluon emission in relativistic plasmas, JHEP 06 (2002) 030 [hep-ph/0204343] [INSPIRE].

[5] G.D. Moore and N. Schlusser, The nonperturbative contribution to asymptotic masses, Phys. Rev. D 102 (2020) 094512 [arXiv:2009.06614] [INSPIRE].

[6] J. Ghiglieri, G.D. Moore, P. Schicho and N. Schlusser, The force-force-correlator in hot QCD perturbatively and from the lattice, JHEP 02 (2022) 058 [arXiv:2112.01407] [INSPIRE].

[7] S. Caron-Huot, On supersymmetry at finite temperature, Phys. Rev. D 79 (2009) 125002 [arXiv:0808.0155] [INSPIRE].

[8] T. Gorda et al., Soft photon propagation in a hot and dense medium to next-to-leading order, Phys. Rev. D 107 (2023) 036012 [arXiv:2204.11279] [INSPIRE].

[9] A. Ekstedt, Two-loop hard thermal loops for vector bosons in general models, JHEP 06 (2023) 135 [arXiv:2302.04894] [INSPIRE].

[10] A. Ekstedt, Propagation of gauge fields in hot and dense plasmas at higher orders, arXiv:2304.09255 [INSPIRE].
[11] T. Gorda, R. Paatelainen, S. Säppi and K. Seppänen, *Soft gluon self-energy at finite temperature and density: hard NLO corrections in general covariant gauge*, JHEP 08 (2023) 021, [arXiv:2304.09187] [INSPIRE].

[12] E. Braaten and R.D. Pisarski, *Simple effective Lagrangian for hard thermal loops*, Phys. Rev. D 45 (1992) R1827 [INSPIRE].

[13] A.D. Linde, *Infrared Problem in Thermodynamics of the Yang-Mills Gas*, Phys. Lett. B 96 (1980) 289 [INSPIRE].

[14] H.A. Weldon, *Covariant Calculations at Finite Temperature: The Relativistic Plasma*, Phys. Rev. D 26 (1982) 1394 [INSPIRE].

[15] S. Caron-Huot, *O(g) plasma effects in jet quenching*, Phys. Rev. D 79 (2009) 065039, [arXiv:0811.1603] [INSPIRE].

[16] E. Braaten, *Solution to the perturbative infrared catastrophe of hot gauge theories*, Phys. Rev. Lett. 74 (1995) 2164 [hep-ph/9409434] [INSPIRE].

[17] E. Braaten and A. Nieto, *Effective field theory approach to high temperature thermodynamics*, Phys. Rev. D 51 (1995) 6990 [hep-ph/9501375] [INSPIRE].

[18] E. Braaten and A. Nieto, *Free energy of QCD at high temperature*, Phys. Rev. D 53 (1996) 3421 [hep-ph/9510408] [INSPIRE].

[19] K. Kajantie, M. Laine, K. Rummukainen and M.E. Shaposhnikov, *Generic rules for high temperature dimensional reduction and their application to the standard model*, Nucl. Phys. B 458 (1996) 90 [hep-ph/9508379] [INSPIRE].

[20] K. Kajantie, M. Laine, K. Rummukainen and M.E. Shaposhnikov, *3-D SU(N) + adjoint Higgs theory and finite temperature QCD*, Nucl. Phys. B 503 (1997) 357 [hep-ph/9704416] [INSPIRE].

[21] G.D. Moore and N. Schlusser, *Transverse momentum broadening from the lattice*, Phys. Rev. D 101 (2020) 014505 [Erratum ibid. 101 (2020) 059903] [arXiv:1911.13127] [INSPIRE].

[22] G.D. Moore, S. Schlichting, N. Schlusser and I. Soudi, *Non-perturbative determination of collisional broadening and medium induced radiation in QCD plasmas*, JHEP 10 (2021) 059 [arXiv:2105.01679] [INSPIRE].

[23] S. Schlichting and I. Soudi, *Splitting rates in QCD plasmas from a nonperturbative determination of the momentum broadening kernel $C(q \perp)$*, Phys. Rev. D 105 (2022) 076002 [arXiv:2111.13731] [INSPIRE].

[24] J. Ghiglieri et al., *Hard parton dispersion in the quark-gluon plasma, non-perturbatively*, PoS HardProbes2023 (2023) 149 [arXiv:2307.09297] [INSPIRE].

[25] E. Weitz, *Theoretical Developments for Jets in Heavy-Ion Collisions*, Ph.D. thesis, SUBATECH, Université de Nantes, 44307 Nantes, France (2023) [arXiv:2311.04988] [INSPIRE].

[26] V.V. Klimov, *Spectrum of Elementary Fermi Excitations in Quark Gluon Plasma* (in Russian), Sov. J. Nucl. Phys. 33 (1981) 934 [INSPIRE].

[27] V.V. Klimov, *Collective Excitations in a Hot Quark Gluon Plasma*, Sov. Phys. JETP 55 (1982) 199 [INSPIRE].

[28] H.A. Weldon, *Effective Fermion Masses of Order $gT$ in High Temperature Gauge Theories with Exact Chiral Invariance*, Phys. Rev. D 26 (1982) 2789 [INSPIRE].

[29] S. Caron-Huot and G.D. Moore, *Heavy quark diffusion in QCD and $N=4$ SYM at next-to-leading order*, JHEP 02 (2008) 081 [arXiv:0801.2173] [INSPIRE].
[30] T. Liou, A.H. Mueller and B. Wu, *Radiative $p_⊥$-broadening of high-energy quarks and gluons in QCD matter*, Nucl. Phys. A 916 (2013) 102 [arXiv:1304.7677] [INSPIRE].

[31] J.-P. Blaizot, F. Dominguez, E. Iancu and Y. Mehtar-Tani, *Probabilistic picture for medium-induced jet evolution*, JHEP 06 (2014) 075 [arXiv:1311.5823] [INSPIRE].

[32] J. Ghiglieri and E. Weitz, *Classical vs quantum corrections to jet broadening in a weakly-coupled Quark-Gluon Plasma*, JHEP 11 (2022) 068 [arXiv:2207.08842] [INSPIRE].

[33] J. Ghiglieri et al., *Next-to-leading order thermal photon production in a weakly coupled quark-gluon plasma*, JHEP 05 (2013) 010 [arXiv:1302.5970] [INSPIRE].

[34] J. Ghiglieri, G.D. Moore and D. Teaney, *QCD Shear Viscosity at (almost) NLO*, JHEP 03 (2018) 179 [arXiv:1802.09535] [INSPIRE].

[35] J. Ghiglieri, G.D. Moore and D. Teaney, *Second-order Hydrodynamics in Next-to-Leading-Order QCD*, Phys. Rev. Lett. 121 (2018) 052302 [arXiv:1805.02663] [INSPIRE].

[36] M. Laine and A. Vuorinen, *Basics of Thermal Field Theory*, Springer (2016) [DOI:10.1007/978-3-319-31933-9] [INSPIRE].

[37] J.I. Kapusta and C. Gale, *Finite-temperature field theory: Principles and applications*, Cambridge University Press (2011) [DOI:10.1017/CBO9780511535130] [INSPIRE].

[38] M. Nishimura and Y. Schröder, *IBP methods at finite temperature*, JHEP 09 (2012) 051 [arXiv:1207.4042] [INSPIRE].

[39] A.I. Davydychev and Y. Schröder, *Recursion-free solution for two-loop vacuum integrals with “collinear” masses*, JHEP 12 (2022) 047 [arXiv:2210.10593] [INSPIRE].

[40] P.B. Arnold and C.-X. Zhai, *The three loop free energy for high temperature QED and QCD with fermions*, Phys. Rev. D 51 (1995) 1906 [hep-ph/9410360] [INSPIRE].

[41] N. Brambilla, J. Ghiglieri, P. Petreczky and A. Vairo, *The Polyakov loop and correlator of Polyakov loops at next-to-next-to-leading order*, Phys. Rev. D 82 (2010) 074019 [arXiv:1007.5172] [INSPIRE].

[42] P.B. Arnold and C.-X. Zhai, *The three loop free energy for pure gauge QCD*, Phys. Rev. D 50 (1994) 7603 [hep-ph/9408276] [INSPIRE].

[43] J. Moller and Y. Schröder, *Open problems in hot QCD*, Nucl. Phys. B Proc. Suppl. 205-206 (2010) 218 [arXiv:1007.1223] [INSPIRE].

[44] J. Ghiglieri and D. Teaney, *Parton energy loss and momentum broadening at NLO in high temperature QCD plasmas*, Int. J. Mod. Phys. E 24 (2015) 1530013 [arXiv:1502.03730] [INSPIRE].

[45] J. Ghiglieri, A. Kurkela, M. Strickland and A. Vuorinen, *Perturbative Thermal QCD: Formalism and Applications*, Phys. Rept. 880 (2020) 1 [arXiv:2002.10188] [INSPIRE].

[46] M. Laine and Y. Schröder, *Two-loop QCD gauge coupling at high temperatures*, JHEP 03 (2005) 067 [hep-ph/0503061] [INSPIRE].