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

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.


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 lowerenergy 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][2][3] for recent reviews.These selfcollimated 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 ⊥ ) .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 ∞ , the asymptotic mass.
The transverse scattering kernel C(b ⊥ ) and m ∞ 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 ∞ , 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] where m 2 ∞,q applies for quarks and m 2 ∞,g applies for gluons.Here C F = (N 2 c − 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 = 1 2 .Eq. (1.1) can be understood as arising from integrating out the energy scale of the jet E ≫ 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 T rather than E ≫ 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][9][10][11], namely the two-loop and power corrections to Hard Thermal Loops and their relation to the asymptotic mass for quasiparticles obeying T ≫ p ≫ 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] where v µ = (1, 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 . . .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 ∼ gT , and gluonic ultrasoft modes, with momenta of O(g 2 T ).The former cause the perturbative expansion to be in g rather than α 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) ≈ T /p ≫ 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.
where m 2 D is the leading-order Debye mass, given by [14] m 2 D = g 2 T 2 3 (C A + T F N F ) + O(g 4 ) . (1.5) The next-to-leading order (NLO) correction, δZ g , is negative and, for values of the coupling corresponding to QGP temperatures, of comparable magnitude to the leading-order (LO) contribution Z LO g .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 α 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][17][18][19][20], a dimensionallyreduced Effective Field Theory (EFT) of QCD.This allows for non-perturbative evaluations of the EQCD contribution to C(b ⊥ ) and Z g through lattice EQCD.We refer to [21][22][23] for the evaluation of C(b ⊥ ) in EQCD, its merging with the T -scale contribution and its impact on medium-induced radiation.
For m ∞ , 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 LO g .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 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 Sec. 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 Sec. 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 Sec. 5, we summarise our main findings and draw our conclusions.The appendices are dedicated to the conventions and technical details of our calculations.

Background
As we anticipated, Z g can be rewritten as the following integral over light-like separated fieldstrength insertions on a light-like Wilson line where the Wilson line itself reads U R (x + ; 0) = P exp ig after having fixed v µ = (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 Fig. 1.It yields [26][27][28] where n B (q 0 ) ≡ (exp(q 0 /T ) − 1) −1 is the Bose-Einstein distribution.Dimensional regularisation in D = 4 − 2ǫ (d = 3 − 2ǫ for the spatial dimensions) has been used to remove the scale-free vacuum contribution proportional to θ(q 0 ).See Appendix A and Eq.(A.2) for our conventions.

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) ≈ T /q, so that Eq. ( 2.3) has an order-g sensitivity to the soft scale gT , viz.q∼gT T q 2 ∼ 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] scale T scale gT scale + O(g 3 ) . (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 O(g 0 ) from Eq. (2.3) is displayed in (2.4) and the first O(g) correction from soft modes in (2.5).At the next order, all scales contribute: the completely non-perturbative color-magnetic scale g 2 T , as well as the soft modes of the color-electric scale gT , and the thermal modes from T .The coefficients c ln hard and c ln soft 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 min : for the moment we have kept the notation in Eqs.(2.4)-(2.6)generic, with µ h and µ 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-depedent.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 ∞ at 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 min dependence, to dimensional regularisation.We will also compute the T -scale contribution in Eq. (2.6), c ln hard 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].

EQCD setup
One starts from the EQCD continuum action, where Φ 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) (2.9) Since the field operators are classical in EQCD, they commute and allow us to write (2.10) Thus, the EQCD equivalent of Z g reads, given rotational invariance in the transverse plane where we used At LO in perturbative EQCD, we can connect the operators with EQCD propagators to find [7] where the massive propagator is the Φ 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) q 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 µ h in Eqs.(2.4) and (2.5).Thus, dimensional regularisation corresponds to subtracting off that (d − 1) q 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 which leaves us with where we have implicitly defined F F ≡ EE − BB − i EB .
In [5,6] continuum-extrapolated lattice EQCD determinations of the three components of F F have been performed for distances 0.25 < g 2 3d 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 ≪ 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 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 µ 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 3d g , [6] merged together the perturbative, intermediate and IR regions as follows where [5,6] The shortest and longest separations for which lattice data is available are L min = 0.25/g 2 3d and L max = 3.0/g 2 3d . 5Between these two separations, the lattice-determined F F lat is used.At asymptotically large distances the modeled IR tail F F tail is used and at short distances the NLO perturbative EQCD result F F 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 Sec. 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 .

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 Fig. 1.As we are dealing with 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 A − A − or A − A x , both of which vanish in Feynman gauge.Furthermore, following an argument to be given in Sec.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.

Equations of motion of the light-like Wilson line
Looking at Eq. (2.1), we can write where Continuing onwards, let us recast Eq. (3.1) as a double integral, using translation invariance where we have shifted x + in going to the second line.We now rewrite Z g as with the index labelling to be elaborated on momentarily.By calling on a method introduced in [29], we rewrite 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.We start with the first term from Eq. (3.3), For the second term, we have, employing Eq. (3.5) Translation invariance dictates that for the first part, so that, by aptly redefining the integration variable, one finds, including its counterpart, that where we used the fact that, in a non-singular gauge, A ⊥ (∞) = 0 and we applied translation invariance and relabeled y + into x + .Finally, for the term with both commutators, we find 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 ⊥ 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 ⊥ ones.However, that would give rise to the following colour structure which vanishes, due to the symmetric (anti-symmetric) nature of A ⊥ a A ⊥ 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].

Decomposition of the gluon self-energy
We thus dive straight into the computation of diagram (c) 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 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 Its retarded and advanced components read 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 We can then rewrite Eq. (3.15) as which leaves us with the decomposition where adv.stands for the advanced part while 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.
For the computation of Z (c) g,q −2 , we can resort to standard Euclidean techniques.We start from where we have kept the even-in-q 0 part of the integrand only.By analyticity, we have Here, Q 2 E ≡ q 2 0 + q 2 is the Euclidean four-momentum squared, with q 0 = 2πnT a bosonic Matsubara frequency and Σ Q E the standard sum-integral, see Eq. (A.2).Likewise, Π E T (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)) where B stands for bosonic, γ E is the Euler-Mascheroni constant and A ≈ 1.28243 is the Glaisher constant.Here and in all following formulae we truncate the expansion in D = 4 − 2ǫ to the order-ǫ 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 Sec. 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)  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 β 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.
Let us then go back to the leading-order term in Eq. (2.3), recalling that, as per Eq.(1.1), m 2 ∞ is obtained by multiplying Z g by g 2 . 6Hence, if that g 2 is replaced by Eq. (3.26), the O(g 4 ) term there generates the counterterm  Here, we relabeled μ → μUV to distinguish it from other occurrences of μ, which are to be understood as factorisation scales, whereas μ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.

Z (c) g,L-T 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ε) 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ǫ 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 Sec. 4 for the L − T part of the soft contribution.
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 selfenergies 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.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 While the vacuum term would be scale-free and vanishing in DR, it physically corresponds to a UV divergence for q z → ∞, which is accompanied by a collinear divergence.This shows up for 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 ⊥ and large q z regions.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]. 9The large-q z region is instead likely related to the emergence of the q z ∼ E region, with E ≫ 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. θ(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 where 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 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 ⊥ /(2k)) must come from two regions of the q z integral in Eq. (3.34), once we take the asymptotic UV expansion q ⊥ ≫ k.Indeed, we have a first region for q z ∼ q ⊥ ≫ k and a second region for q z ∼ q 2 ⊥ /k ≫ q ⊥ -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 and separating them which reproduces the q ⊥ ≫ 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 + ∼ E region.Its proper handling is also left to future work.We conclude by remarking that in the real-time calculation of Z (c) g,q −2 (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 (c) g,q −2 -see [25] -might explain why no such cancellation happens for Z (c) g,L-T .

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.
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/ǫ 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 Sec.3.1 we then know that this c 2 -proportional UV divergence must come from the EQCD analogue of diagram (c), shown in Fig. 2.This suggests that where Z 3d (c) g,q −2 ,c 2 is the (q − ) 2 and Z 3d (c) g,L-T ,c 2 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 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 ) where δZ non-pert class+NLO T g 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 3d (c) g,q −2 ,c 2 and Z 3d (c) g L-T ,c 2 , starting from the former.With respect to the decomposition defined in Eq. (2.16), we first take the BB contribution from Appendix A.3 of [6] (cf.Eq. (A.14)) where Π T A a i A b j div is the EQCD transverse gluon polarisation tensor in the massless limit.As we are investigating a UV divergence on the EQCD side, the q ≫ m D limit is the appropriate one.It can be extracted either from the m D → 0 limit of Eq. (A.10) in [6] or from the zero-mode contribution to Eq. (B.4).We find where we have denoted the d = 3 limit for future use.
The second term in Eq. (4.7), proportional to (d − 1)q 2 z , corresponds to the Z 3d (c) g,q −2 ,c 2 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 match g contribution undoes the corresponding piece of the c 2 subtraction.Looking at Eq. (2.16) and using Eq.(4.8), we then find Since g 2 3d = g 2 T + O(g 4 ), the 1/ǫ term cancels against the one from Eq. (4.1).As for the L − T piece, we need the q 2 ⊥ -proportional piece of Eq. (4.7) together with the contribution from the EE correlator.We extract this from Eq. (A.8) of [6] and expand towards the UV to arrive at As the corresponding IR pole has been obtained using the mixed scheme described at the beginning of Sec.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 To get the full Π L − Π T EQCD contribution, we need to add to this its BB counterpart.It arises from the q 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) which, reassuringly contains a 1/ǫ 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 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 EE and a BB one, to check the explicit cancellation of 1/ǫ 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 2 3d = g 2 T , to find Note that in computing Z non-pert class g , the ln L min above cancels against the − dL L c 2 g 2 3d 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 3d g 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 Sec. 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 The numerical impact of Eqs.(4.14) and (4.15) will be discussed in the next section.

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 Sec. 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 (5.1) The one outstanding roadblock is that our evaluation of the "longitudinal minus transverse" Z (c) g,L-T n =0 contribution in Eq. (3.29) contains yet uncancelled double-logarithmic divergences from the sensitivity to the collinear (q ⊥ ∼ m ∞ ) and hard (q + ∼ E ≫ 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 (c),IR/UV safe g,q −2 , 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 which is obtained by summing Eqs.(3.28) and (4.3).We also set μUV to Eq. (D.1), in keeping with the renormalisation scale used in the EQCD matching, and L min = 1/(4g 2 3d ), the value used in Z 3d g 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 (c) IR/UV safe g,q −2 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 (c) g,L-T ,n =0 , 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 3d g merge from [6].δZ finite g Eq.(4.15) is defined by replacing Eq. (4.14) with Eq. (4.15) in Z non-pert class g , 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 g,LO 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 LO g = T 2 /6, as per Eq.(4.5).
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 δZ finite g are negative and, in magnitude, larger than Z LO g = 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 3d g,LO in the last column of the table.Recall that Z 3d g,LO 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 δZ finite g,Eq.(4.15) 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 + ∼ E ≫ 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 Sec.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 ≫ T when the T -scale contribution to these operators is considered.If the soft contri-bution is considered, the expansion becomes E ≫ m D , so that, in the weak-coupling regime, it remains applicable for E ∼ 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. 10This has then the very important consequence that our values for Z non-pert class g , tabulated in Tab. 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][34][35].

A. Conventions
Our sign for the covariant derivative is 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 where the two light-like reference vectors are defined as vµ ≡ 1 2 (1, 0, 0, −1) , v µ ≡ (1, 0, 0, 1) . 10We remark that, even in this regime, our results remain different from those of [8][9][10][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 ≪ Qi ∼ T , where P is the external momentum and Qi the loop momenta, to then take the P ≫ mD limit.Our calculation, on the other hand, corresponds to taking a P ≫ Qi ∼ T expansion to then take the P → T ≫ mD 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][9][10][11] are finite after charge renormalization, further highlighting their difference.
This asymmetric convention for the + and − components of the light-cone coordinates has two advantages: it has unitary Jacobian, i.e. dp 0 dp z = dp + dp − , in scalings where p − ≪ p + it then implies p 0 ≈ p z ≈ p + .

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 Sec. 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) The main ingredient for the Euclidean evaluation of the (q − ) 2 contribution is the Feynmangauge transverse Euclidean polarisation tensor, which can be obtained from, e.g.[36].It reads 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.
where B and F stand for bosonic and fermionic.

C. Diagram (c) computation: technical details
We start from Eq. (3.22).Upon inserting Eq. (B.4) into Eq.(3.22), we encounter the following class of sum-integrals where I 00 0;s 1 s 2 s 3 = I s 1 s 2 s 3 and fermionic momenta are denoted by curly brackets.The integrals that appear in Eq. (B.4) are and the fermionic counterparts Ĵi obtained by summing over odd {K E }.We used the one-loop bosonic and fermionic master integrals where ζ s = ζ(s) is the Riemann zeta function.J 2 vanishes identically, as shown with integrationby-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.
The fermionic counterparts read We used the one-loop fermionic master integrals (C.9) and that Ĵ2 too vanishes identically, see e.g.[40].Details on Ĵ3 and Ĵ5 are provided below in Sec.C.1.1.
From Eq. (3.22), Eq. (B.4) and the definitions of the master integrals, we have The bosonic part gives 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) C.1.1.Evaluation of J 3 and J 5 and their fermionic counterparts For the bosonic J 3 and fermionic Ĵ3 , we first perform the Matsubara sums and drop the scale-free vacuum × vacuum term, i.e.
where we have used the q ↔ k symmetry of the integrand.We can now split the vacuum and thermal parts in k.The former can be treated as follows 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 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.
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 k integration followed by the q one.We find Using the position-space method of [40,42], the numerical part can be integrated analytically with full agreement, finding To evaluate J 5 , we employ again the positon-space methods of [42] and write J 5 as where we can separate the vacuum and thermal contribution to Π(Q) The vacuum contribution is obtainable by replacing the Matsubara sum over k 0 with a Euclidean dk 0 /(2π).Employing standard DR methods and inserting the result into Eq.(C.32), yields (C.34)For Π (T ) , we instead proceed in position space, i.e.
where δ p 0 = δ p 0 ,0 and in the last step we took the three-dimensional Fourier transforms.One can also work out the Matsubara sum, finding where q0 ≡ q 0 /(2πT ) and r ≡ 2πT r.To get Π (T ) , we must subtract the vacuum contribution, again obtained by replacing the sum with the integral.We then find In principle, we should naively insert this into Eq.(C.32).We would find that 1) Since the Fourier transform 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 ≫ T limit of Eq. (C.37) in position space and add it back in momentum space in DR.
We then find where the dots stand for higher-order terms in the expansion for r ≪ 1 with r q0 ∼ 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 By splitting J (T ) 5 into three terms and noting that primed sum-integrals exclude the zero mode, we obtain For the first term we proceed in position space.After using Eq.(C.38) and performing the q 0 sum, we have where the integration is finite and has been carried out numerically. 11 For the divergent piece we find, using IBP 11 An analytical evaluation could be possible by using IBP: ln(1 + e −r ) = ln(1 − e −2r ) − ln(1 − e −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 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 where s {i} = j∈{i} s j .The two special cases relevant for our computation

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 which, after inserting Eq. (C.52) into (C.51),gives rise to a pole and a cut contribution.The pole contribution is where in DR, the second step used that the vacuum part is scale-free.The cut contribution of the HTL piece is instead 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 (C.56) Our strategy here is to do the following subtraction: n B (q 0 ) → n B (q 0 ) ∓ (T /q 0 − θ(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 For the DR subtracted part, we instead do the frequency integral first, splitting it into the q > T and q < T ranges, where the principal value (PV) P arises from the ε → 0 limit.The q > T slice is finite whereas the T > q part is IR-divergent and needs regularisation Let us now consider the non-HTL piece where δΠ = Π − Π HTL .This then gives Here, we used the fact that the 1 2 piece vanishes 12 to restrict to the purely odd 1 2 + n B (q 0 ) component.Let us again symmetrise and restrict to twice the positive frequency range We start by dealing with the vacuum part 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 ε → 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 ε → 0 limit.We can thus proceed to the q z integral in the two ranges 2k > q and q > 2k.

(C.70)
As was done in Sec.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ε) 2 factor, the q 0 integrand possesses poles in addition to those associated with the Matsubara frequencies, meaning that we cannot effectively trade the q 0 integral for a sum over Matsubara modes as was done before.However, by first shifting: q z → q z + x 0 x z q 0 , we may eliminate the q 0 pole tied to the 1/(q − − iε) 2 prefactor, before proceeding to carry out the q 0 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] 13 that as long as X is lightlike or spacelike, one may always work in a frame where |x 0 | |x z | ≤ 1, in which case the analytical structure of the causal functions is not altered. 14Thus, in this case where X is lightlike, the appropriate shift is precisely: q z → q z + q 0 .Following this maneuver, we can safely carry out the integral q 0 and extract the n = 0 contribution which is needed for the matching calculation completed in Sec. 4.

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 Tab. 2, together with the data for Z 3d g merge from [6]. 16The matching procedure corresponds to a renormalisation scale μUV = 4πT e −γ E μ , (D.1) where at one-loop level [20] μ = exp Hence, we set μUV to Eq. (D.1) when numerically evaluating Eq. (3.28) in Sec. 5.

Figure 1 :
Figure 1: Diagrams contributing the QCD force-force correlator Z g (2.1) at leading and next-toleading order.An external gray shaded vertex denotes an F −⊥ insertion; internal 2-point blobs the respective self-energy; a solid line a Wilson line; and a curly line a gauge boson (A µ ).

. 5 )
This should shed light on the labelling in Eq. (3.3): the indices there denote whether the −∂ ⊥ A − ( ⊥ ) or the [D − , A ⊥ ] component ( − ) of the field strength tensor has been taken.

Figure 2 :
Figure 2: EQCD equivalent of diagram (c).The labelling is identical to that from Fig. 1 except that here, curly lines represent spatial gauge bosons, double lines are adjoint Wilson lines and solid lines are adjoint scalars.
T n =0 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 (c) g L-T n =0 , 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 2 −
-T ,F,non-HTL th =g 2 T F N F T 2 48π -T ,F = g 2 T F N F T 2 48π 2 make the replacement 2N F T F → C A in Eq. (C.61) to obtain the bosonic equivalent.This yields Z

For 2
the vacuum part of the non-HTL contribution, we have instead Z (c) g,L-T ,B,non-HTL vac = g

Z 2 C
-T ,B,non-HTL th n =0 = g Q zero-mode contribution, the total is thenZ (c) g,L-T ,B,n =0 = g 2 C A T 2 48π 2 = ζ(s) is the Riemann zeta function, (ln ζ s ) ′ = ζ ′ (s)/ζ(s) and γ 1 is the first Stieltjes constant.For the C A -proportional gauge contribution, we have again subtracted off the zeromode 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/ǫ 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)

Table 2 :
No zero-mode subtraction was necessary in the HTL part, since the zero mode is scale-free there.