The force-force-correlator in hot QCD perturbatively and from the lattice

High-energy particles traversing a medium experience modified dispersion. In the Quark-Gluon Plasma, such dispersion affects jet propagation and transport properties and should be determined better. Above $\sim 2T_c$ we expect strongly coupled infrared behavior and perturbative ultraviolet behavior, allowing a perturbative matching to an effective theory called EQCD, which can be studied non-perturbatively. We study the relevant non-local operator in EQCD at next-to-leading order which allows for a complete EQCD-to-lattice match and prepares the groundwork for a matching between EQCD and full QCD. Our results in EQCD show remarkable agreement between perturbation theory and the lattice in the expected regime.


Introduction
The Quark-Gluon Plasma (QGP), an exotic state of strongly interacting matter, is currently investigated in heavy-ion collision experiments. Its experimental characterization and the accompanying theoretical activity proceed along the axes of bulk properties and hard probes. The former is the study of the behavior of the many lower-energy produced particles, which are understood to arise from the hadronization of a hydrodynamically-evolving, near-equilibrium medium, while the latter concentrates on the few particles which are very energetic or weakly coupled to the medium.
Jets, a key hard probe, are important experimental sources of evidence about the nature of the strong nuclear interaction under extreme conditions -see [1,2] for recent reviews. Jets are generated by colored particles at high energies, a regime where the theory of the strong nuclear force, Quantum Chromodynamics (QCD), is supposed to be weakly coupled and accessible to perturbative methods. Here we focus on jets created by light constituents such as light quarks (u, d, s) or gluons. It was found by Klimov [3,4] and independently by Weldon [5,6] that -though being massless -high-energy particles with momentum k follow the dispersion relation of massive particles when traversing the QGP (1.1) By scattering with the medium they acquire an effective mass, which, at very large momentum, is called the asymptotic mass. To one-loop order these masses are composed of the gauge condensate Z g and the fermion condensate Z f : 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 is the Dynkin index for the quarks. The condensates Z g and Z f are non-local and have a gauge-invariant definition 1 in terms of correlators [7,8] where v µ = (1, v) is the light-like four-velocity of the hard particle, d R,A are the dimensions of the fermion and adjoint representations respectively, and the expectation value . . . denotes a thermal expectation value. Our conventions are that the metric is the mostly-plus one, the covariant derivative is D µ = ∂ µ − igA µ (with the gauge coupling g) and the field strength tensor is F µν = i g [D µ , D ν ]. Since we are interested in QCD, henceforth we will specialize to the gauge group SU(3), with C F = 4/3, C A = 3, d R = 3, d A = 8, and T F = 1/2 where not indicated differently.
Eq. (1.4) can be rewritten in coordinate space, where the inverse powers of derivatives describe an integral over separations, with a Wilson line to reflect that they are covariant derivatives. That is, Z g is an integral over the lightlike separation x + ≡ (x 0 + v · x)/2 of a correlator of two covariant Lorentz-force insertions (1.5) where U A (x + ; 0) is an adjoint, lightlike Wilson line. 2 Despite the jet momentum being much larger than the temperature of the medium that it traverses, the interaction between jet and medium still receives contributions from the infrared (IR) regime, i.e. from energy-momentum regions of O(gT ) or smaller. These contributions can, for many quantities such as the interaction rate, be dominant. Indeed, the emergence of these IR scales causes the perturbative series to be an expansion in g rather than in α s ; furthermore, at an operator-dependent order in g, the perturbative expansion breaks down, due to the presence of the so-called magnetic or non-perturbative scale g 2 T [9]. In our case, the IR contributions affect the two operators at different orders. The fermionic condensate Z f is well under perturbative control, with no IR contribution appearing as an O(g) correction to the leading order (LO) result [8] (1.6) Only at higher orders it receives contributions from the IR regime. However, the gauge condensate Z g receives IR contributions already at O(g): the next-to-leading order (NLO) result reads [8] Z g = Z LO g + δZ g = where at leading order the Debye screening mass m 2 D is given by [5] (1.8) Formally m 2 D is a parameter of an infrared effective description, EQCD, which we will introduce momentarily. At leading order it equals twice the gluonic asymptotic mass, m 2 D = 2m 2 ∞,g . For the values of T and therefore g which are relevant in any conceivable heavy-ion experiment, the Debye mass is large enough for the two known terms in (1.7) to be comparable. Therefore, the IR contribution spoils the convergence of the perturbative expansion. 2 The Wilson line structure of (1.5) is a somewhat delicate issue that requires further clarification. Grouptheoretically, an adjoint line is equivalent to a pair of fundamental Wilson lines with the same endpoints, but they may differ in terms of operator ordering or endpoints, as the fundamental lines may stretch back to x + = −∞. If we treat v = 1 + ǫ then all operators along the Wilson line commute, and there is no difference between the treatments. But if collinear physics becomes important at some higher order in the 4-dimensional treatment, then there will be a discontinuity between a v = 1 + ǫ and a v = 1 − ǫ treatment and the distinction between Wilson line types may become important. This cannot occur in eq. (3.4), the EQCD version of eq. (1.5) which we study in this paper, because of the time-independence of EQCD. We will come back to this issue when addressing the T -scale contribution at O(g 2 ) in an upcoming paper.
In a remarkable set of papers [8,10] it was however shown that the leading soft interactions of a very hard, hence light-like, parton with a soft thermal bath can be greatly simplified, both for the interaction rate [10] and for Z g [8]. In more detail, they can be isolated and treated in a dimensionally-reduced Effective Field Theory of thermal QCD, called Electrostatic QCD (EQCD) [11][12][13][14][15]. Its effective parameters encode the ultraviolet (UV) behavior of thermal QCD as a function of temperature T and number of massless quark flavors N f . This is achieved by matching the Green's functions of both theories in the IR. To investigate a QCD operator using this effective field theory, we must determine the EQCD counterpart of that operator; for us, this means determining the EQCD equivalent of the correlator of eq. (1.5).
Our long-term goal, started in [16], is to determine Z g by a matching between full QCD and EQCD and by computing the EQCD contribution non-perturbatively on the lattice. This paper will take two steps in this direction, leaving one step for future work. First, we will update the lattice part of the EQCD equivalent to eq. (1.5). Second, we will calculate the EQCD correlator to next-to-leading order in EQCD perturbation theory. This improves the fitting necessary in the lattice determination, and it also sets the groundwork for the matching to full QCD, by gaining a better analytical understanding of the (unphysical) UV divergences in EQCD, and how these divergences are appropriately dealt with. We argue on dimensional grounds that NLO is the last order at which potentially UV-divergent terms from EQCD in (1.5) can arise. When supplied with a full four-dimensional perturbative calculation at the same order, the matching of our lattice EQCD results back to full QCD would be complete, and entirely IR-resummed values for m 2 ∞ can be provided. The four-dimensional calculation is spared for a separate publication, however. Moreover, one can finally only interpret the UV divergences in EQCD in the context of the full four-dimensional framework. Therefore, a fully consistent subtraction scheme also has to be postponed to another publication. This paper is organized as follows: Section 2 clarifies the power counting scheme(s) used throughout this present publication. Section 3 details the transition from full QCD to EQCD. We present the analytical calculation of the force-force correlator in EQCD at NLO in sec. 4 and update the existing lattice calculation in sec. 5. Section 6 discusses our results with an outlook to future interesting research. Details of our perturbative calculation are found in appendix A.

Power-counting scheme
To shed light on the at first sight opaque power counting scheme of our result, we find it instructive to parametrize our result before actually diving into the computation. Let us start with a power-counting analysis of Z g . Its leading order form is where n B is the Bose-Einstein distribution. This form shows how the LO term in eq. (1.7) is obtained: it is the leading contribution from the momentum region p ∼ T . The first IR contribution arises from p ∼ gT , which is Bose-enhanced, i.e. n B (p) ≈ T /p ∼ 1/g. This, together with d 3 p/p ∼ g 2 T 2 , suggests an O(g) contribution from this scale, corresponding to the contribution to Z g from the leading order in EQCD. Going further to the IR, we can expect the first contribution from the non-perturbative scale g 2 T at O(g 2 ). This is precisely the order reached by our NLO-in-EQCD perturbative calculation, as well as that of the second-order contribution from the scale p ∼ T . Hence, a strict perturbative expansion of Z g in the QCD coupling g would yield where the three lines correspond to O(g 0 ) -O(g 2 ), respectively, whereas columns correspond to the originating scale. At O(g 2 ) all scales contribute, and the contribution of each of them is scheme-dependent: we indicate this through the intermediate regulators T ≫ µ h ≫ gT and gT ≫ µ s ≫ g 2 T . 3 These should not be taken as indicating the choice of a cutoff scheme, but rather as generic placeholders for an arbitrary choice of scheme to separate the contribution of the scales. This scheme dependence also affects the coefficients c T , c gT and c g 2 T , with the latter being non-perturbative. Through our calculation we will be able to determine c ln hard and c ln soft from the logarithmic divergences of the NLO-in-EQCD calculation. 4 That is not, however, our main aim, which is instead to incorporate the all-order EQCD contribution to Z g through lattice EQCD.
Finally, we remark that, due to the super-renormalizable nature of EQCD, the EQCD contributions of O(g 3 ) and higher will not present further UV divergences. In the language of this section, we can expect the UV behavior at O(g n ) to be µ 2−n h , so that, after the leading linear and subleading log divergences we can only find negative power laws in the UV cutoff.
3 Light-cone observables and dimensionally-reduced theories As we noted, it was found by Caron-Huot that the infrared, non-perturbative part of the jet-medium interaction can be isolated in the framework of electrostatic QCD (EQCD) [10].
Since the gluonic Matsubara zero-mode contributes a factor of 1/g for each closed thermal loop, the perturbative power counting scheme of a g 2 -suppression per loop is spoiled. In other words, the loop and coupling expansions misalign, which is known as Linde's infrared problem [9]. This problem can be bypassed by reorganising the perturbative expansion, which is most economically achieved by treating the zero mode separately in a three-dimensional effective theory [17,18]. The transition from fundamental four-dimensional, thermal QCD to its three-dimensional effective theory EQCD is dubbed dimensional reduction. The continuum action of this effective theory reads 5 The former temporal component of the gauge field turns into the scalar field viz. Φ = iA 0 in the adjoint representation of SU(3) with the acquired thermal mass m 2 D . Since the resulting theory is static in three dimensions, derivatives in temporal direction are absent. Similarly, gauge invariance does not protect the A 0 or Φ field from acquiring a mass. The quartic coupling of Φ is a remnant of the four-A 0 -interaction that arises at twoloop level in full-QCD perturbation theory. This effective theory was originally proposed by Appelquist and Pisarski [17,18]. Braaten and Nieto [12] presented a rigorous perturbative procedure for constructing the effective theory and determining its parameters m 2 D , λ E , and g 2 3d by matching; they also named the theory 'Electrostatic Quantum Chromodynamics' (EQCD). The current matching is available at O(g 6 ) for g 2 3d [19], m 2 D [20] and λ E [15], with further improvements in [21,22]. The matching of these effective parameters and the running coupling g(μ/Λ MS ) 2 are employed at O(g 4 ), while we use a renormalization scale of Λ MS = 341 MeV [23], obtained from (2 + 1)-flavor lattice simulations.
Since EQCD is a super-renormalizable theory, all amplitudes can be rendered finite by a finite number of counterterms. In our case, only the scalar mass m 2 D receives a renormalization and therefore carries the only scale dependence of any parameter in the action. Using the dimensionful gauge coupling to set the scale µ = g 2 3d , one can re-phrase all EFT parameters as dimensionless ratios Concerning the matching of the operators in (1.5), dimensional reduction involves replacing F i0 → i D i , Φ , as mostly already outlined in [16]. Furthermore, we replace the lightlike Wilson line U A with its EQCD counterpart [10] oriented along the z-axis and suppressed constant transverse coordinates (x ⊥ ). With rotational invariance in the transverse plane, we find where the explicit factor of T accounts for the different normalisation of the EQCD fields and the three different correlators are abbreviated as An adjoint Wilson line between two operator insertions O and O ′ can be related to a pair of fundamental Wilson lines via On the lattice it is more convenient to evaluate the correlators contributing to (3.4) in the fundamental representation. By absorbing a factor of 1/2 from eq. (3.4) into the correlators (3.5)-(3.7), they can naturally be considered as living in the fundamental representation. We will approach computing Z 3d g from (3.4) from two different sides: perturbatively to NLO in EQCD or O(g 2 3d ), and non-perturbatively in lattice EQCD. At short separations L, we expect the perturbative and non-perturbative lattice results to agree. Treating short distances on the lattice is challenging. The lattice spacing must be kept several times smaller than the separation of interest to avoid contamination from higher-dimension operators. The lattice volume must be kept larger than a certain physical scale to ensure that one maintains the right phase structure (the symmetric phase of EQCD is only metastable, and the metastability is lost when the volume gets too small [16]). And the precision we demand becomes prohibitive, as the EE and BB correlators each diverge as 1/L 3 , while the errors must be smaller than of order 1/L 2 for eq. (3.4) to converge. At large separations, however, perturbation theory is supposed to become unreliable, and data points directly obtained from lattice simulations or fits of large-g 2 3d L models can provide further insight.
4 Perturbative determination of Z 3d g at NLO in EQCD As motivated above, our task is to match the gluonic EQCD correlator (3.4) to its full QCD counterpart. In fact, the UV contributes to Z g since the latter contains an integral over length Figure 1. Diagrams contributing to leading and next-to-leading order to the EQCD force-force correlator Z 3d g (3.4). An external gray shaded vertex denotes a D x Φ or F xz insertion; internal 2-point blobs the respective self-energy; a double line a Wilson line; and a curly line a spatial gauge boson (A i ). A solid line is a placeholder for either an adjoint scalar (Φ) or a spatial gauge boson (otherwise curly). scales in the domain L ∈ [0, ∞). The ultraviolet (UV) region is precisely where we expect corrections from full QCD that are not included in EQCD. Fortunately, these contributions should be under better perturbative control than the IR regime due to asymptotic freedom. The corresponding diagrams are compiled in fig. 1.
Computing the three correlators in EQCD, EE , i EB , and BB , is only possible at finite values of g 2 3d L on the lattice. Beyond the feasible range of g 2 3d L on the lattice, one has to rely on models. For large g 2 3d L, one can fit the largest-g 2 3d L lattice data points to asymptotic models, as done in [16] and updated in sec. 5. This regime produces a small contribution because the correlators decay exponentially here. For small g 2 3d L, perturbation theory is supposed to work in EQCD. Since the three-dimensional coupling g 2 3d carries mass dimension one, and the correlators carry mass dimension three, dimensional analysis tells us that the tree-level EQCD expressions can go as 1/L 3 in the small-g 2 3d L-limit at worst, whereas the one-loop level can contain g 2 3d /L 2 at worst, and all higher loop levels are O(g 4 3d /L) at worst. The L dL integration leading to Z 3d g in (3.4) can therefore receive UV-divergent contributions only from the O(1/L 3 ) LO terms or from the O(g 2 3d /L 2 ) NLO terms. 6 All higher-order contributions to Z 3d g are short-distance finite. Therefore, a one-loop analytical calculation of the three correlators is not only required quantitatively for increasing the agreement of lattice data and perturbation theory at small g 2 3d L, but also qualitatively for a comprehensive treatment of all possible UV divergences. Note that this short-distance region is where EQCD no longer provides a good description of full QCD. Nevertheless, an accurate treatment of this region will be needed when we carry out the matching to the full four-dimensional theory, which we leave to a future publication.
Below, we present the next-to-leading order perturbative calculation in EQCD. Its mass parameter m 2 D is fully resummed and not treated as a perturbation of O(g 2 ), while also fully taking the quartic Φ-vertex into account. Furthermore, we employ momentum-space gauge 6 When we talk about leading or next-to-leading order in the following, we refer to orders in the EQCD perturbative expansion in g 2 3d , being related but not to be confused with the full QCD perturbative expansion in g. We will thus drop the previously adapted -in-EQCD specifiers.

and adjoint scalar propagators
where the former is in general covariant gauge with gauge fixing parameter ξ and p is a three-dimensional vector with modulus |p| = 0. We will generally use Feynman gauge ξ = 1, and this should be assumed except in expressions with explicit ξ dependence.

The leading order
The leading O(g 0 3d ) contribution is also the free solution to Z 3d g , consists of the graphs from fig. 1(a), and was calculated in [16].
After inserting the propagators from eqs. (4.1) and (4.2), and explicitly separating integration momenta as p = p ⊥ + p z e z , the remaining integrals are Fourier transforms of the form where the gauge dependence of BB cancels due to the antisymmetry of the field strength tensor. The integrals over p z and p ⊥ can be evaluated with the residue theorem and in dimensional regularization respectively with d = 3 − 2ǫ, and results are collected in appendix A.1.
Contracting color indices δ aa = 2C A C F yields for the diagonal correlators with Wilson line endpoints x = (0 ⊥ , L) and x ′ = (0 ⊥ , 0). Here and in the following we draw and evaluate each diagram (e.g. (a)) in its adjoint form which is, as we argued, twice the fundamental form in the definitions of EE , BB , and i EB below eq. (3.5). This simplifies the graphical notation and color evaluation. In turn, we account for the factor of 2 explicitly in contributions of adjoint diagrams to any channel e.g. EE , as 2 × (a) EE . Henceforth, a curly line denotes a spatial gauge boson (A i ) and a solid line an adjoint scalar (Φ). The Z 2 symmetry of EQCD requires the number of Φ fields to be even. 7 Therefore, at leading order the mixed i EB correlator vanishes and at next-to leading order only graphs with a single Φ sourced from the Wilson line contribute. Similarly, only graphs with zero or two Φ fields sourced from the Wilson line contribute at NLO to the EE and BB operators (cf. fig. 1).
The physical short distance behavior of the (four-dimensional) gauge condensate Z g is perturbative, as detailed in sec. 2. Therefore the UV of Z g is described in four-dimensional thermal QCD, as it is dominated by energy scales of order T rather than m D . When it comes to an IR-consistent treatment of longer-range corrections, EQCD comes into play. Since both lattice and perturbative EQCD results for Z 3d g contain an (unphysical) UV limit, while conversely the T -scale contributions contain an incorrect IR limit, as they integrate over the IR momentum region without resummation (see eq. (2.1)), a subtraction scheme is needed to avoid double counting.
The UV limit of the EQCD correlators (3.4), can be extracted from the small m D L asymptote of the perturbative contributions. At leading order starting from (4.4) and (4.5) and taking . . . subtr = lim m D L→0 . . . this gives rise to the subtractions which contain all non-resummed (i.e. m D → 0) EQCD contributions at this order and should be removed from Z 3d g . These subtractions also eliminate the UV divergences in the EQCD result, which are not considered to be physical: they correspond to the linear-in-µ h term in eq. (2.3) (see footnote (3)) and are already included in the unresummed zero-mode contribution to the LO-in-QCD contribution to Z g .

Relation to momentum broadening kernel
We now comment on the relation between Z g and the transverse momentum broadening coefficientq. The latter is defined as the second moment of the broadening probability P (q ⊥ ) or of the transverse scattering kernel C(q ⊥ ) up to some UV cutoff q max . Both can be derived from a Wilson loop in the (x + , x ⊥ ) plane [10,[25][26][27]. Upon writing q 2 ⊥ as transverse derivatives acting on this Wilson loop, one obtains for a fundamental source [27] where U F (b + ; a + ) y ⊥ is a fundamental Wilson line in the + direction from a + to b + at fixed y ⊥ coordinate, while the − coordinate is fixed to the same value for all fields in the expression above.Ū F (b ⊥ ; a ⊥ ) y + is instead a fundamental Wilson line in the ⊥ direction from a ⊥ to b ⊥ at fixed y + coordinate. If the q ⊥ integration is taken with infinite cutoff, i.e. in dimensional regularization, the resulting δ (2) (x ⊥ ) will squeeze the Wilson loop in eq. (4.7) into a Wilson line, viz.
where we dropped the -now fixed -transverse coordinate of the Wilson lines. We can now apply dimensional reduction on eq. (4.8) to find the soft, three-dimensional contribution toq This shows thatq is related to Z 3d g by a change in the overall prefactor and in the integration, which is +∞ −∞ dL rather than +∞ 0 dL L. We stress that eq. (4.9) is to be understood in dimensional regularization only. If this regularization is not taken, divergences arise at L = 0 in the dL integration. This is apparent at LO from eqs. (4.4) and (4.5), which would diverge as dL/L 3 at the origin. We have however checked that, if eqs. (4.4) and (4.5) are evaluated in d dimensions and then plugged in eq. (4.9) and integrated in dL, they yield where our conventions for dimensional regularization are summarized in appendix A.1. This agrees with a dimensionally-regularized evaluation of the LO soft contribution toq [10,28], This evaluation highlights another aspect that will become relevant when performing the same check for our NLO evaluation of Z 3d g , namely that squeezing the Wilson loop onto itself through dimensional regularization obfuscates some cancellations that would be apparent had the integrals been performed in a different order. With that we mean that the UV divergence ofq 3d LO is logarithmic, as it is well known and clear from eq. (4.11). However, inserting eqs. (4.4) and (4.5) in eq. (4.9) leads to severe dL/L 3 UV divergences: the cancellation of the leading UV behavior between the magnetic and electric contributions -made explicit by the two terms in brackets in eq. (4.11) -seems lost. However, let us return to eq. (4.5): had we taken the dL integration before taking the d 3 p integrals in eq. (4.3), the ∂ z ∂ z ′ part would simply vanish due to the δ(p z ) resulting from the dL integration, thereby making apparent the cancellation between the leading UV behavior of EE and BB .

The next-to-leading order
At next-to-leading-order (NLO) multiple graphs contribute to Z 3d g in EQCD as depicted in fig. 1(b)-(j). Relegating their explicit evaluation to appendix A, here we merely sum their individual contributions which yields for the EE -correlator 8 where E 1 is the exponential integral function (A.5) and x = λ E /g 2 3d the dimensionless scalar self-coupling. For the BB -correlator, we find (4.13) Finally, our perturbative prediction for i EB reads where the hyperbolic sine integral function is defined as Shi(x) = g the gauge fixing parameter ξ of eq. (4.1) duly cancels at the integral level.
Another possible source of error is the regularization of divergences. To this end, we extracted UV and IR divergences separately by computing diagrams both in position and momentum space. Where possible we verified individual results by comparison to literature in the soft limit p 2 m 2 D ≪ 1, for instance to [29] for diagram (c). In a final non-trivial crosscheck, we recovered the NLO result for the jet quenching parameterq from [10] by relying on its relation with Z 3d g ; see sec. 4.2. As we discuss there, that relation holds in principle in dimensional regularization only. However, we have not performed the entire NLO determination of Z 3d g in an arbitrary number of dimensions, so we cannot simply take the final dL integral as in eq. (4.9). We have instead worked at the integral level: as the previous subsection discusses, these integrals simplify greatly if the dL integral is taken before the d 3 p ones. Finally, we have been able to exploit the fact that the NLO soft contribution toq, due to the super-renormalizable nature of EQCD, is free of logarithmic divergences. It only contains a linear divergence and is finite in dimensional regularization, so that we have been able to carry out the remaining d 2 p ⊥ integrals and recover the result of [10].
The determination of the UV limit of the EQCD correlators at next-to-leading order changes qualitatively compared to leading order (4.6). The correlator i EB is non-vanishing, and the expressions become more complicated. While the scalar self-coupling λ E enters eq. (4.12) through the scalar self-energy, it fortunately does not contribute any ∼ 1/L 2 terms that would need to be subtracted.
The asymptotic small-m D L limit of the EE EE subtr NLO = 2 indeed diverges in the UV (L → 0) and can therefore be subtracted from the EQCD result for Z 3d g . Correlators BB and i EB are at worst O(L −1 ) in the UV at NLO in EQCD and therefore do not contribute any divergences to Z 3d g and thus require no subtraction. The two correlators containing E-field insertions give an IR-finite contribution to Z 3d g , since they are not directly connected to the (inherently non-perturbative) magnetic sector of QCD at O(g 2 ) of the perturbative expansion, while BB has such direct connection at this order: it behaves as (ln(L) + c)/L 2 for large L and gives thus rise to an IR-divergent contribution to Z 3d g . Lattice EQCD data does not suffer from IR divergences, therefore we are not impacted by perturbative IR divergences as long as we switch from perturbation theory to the lattice at sufficiently small g 2 3d L. However, the necessary subtraction EE subtr NLO in (4.16) additionally introduces, if integrated to infinite L, an IR divergence which will remain necessary for a lattice EQCD treatment due to UV divergences of the lattice data. In the definition of Z 3d g (3.4) the final integration over dL L of O(L −2 ) terms delivers a logarithmic IR divergence which we expect to cancel in the matching with a corresponding IR divergence in full QCD. Consequently, Z 3d g , the result within EQCD presented in this publication, will still depend on an IR cutoff which we will introduce in the next section, even though the dependence on the regulator should vanish in the final, fully matched result for Z g .

EQCD lattice results
In this section, we extend the results of [16] to include i EB and then compare the contiuumextrapolated data to our perturbative determination. We then provide a UV-subtracted, scheme-dependent number for the IR contribution to Z g ; the scheme dependence will disappear once the NLO T -scale contribution will become available. For the details of the used implementation of lattice EQCD we ask the reader to consult appendix B and the references therein.

Lattice determination and continuum extrapolation
Having identified i EB as an additional non-trivial, but numerically subdominant correlation function contributing to Z g , there is a need to measure i EB from the lattice. Just as in the continuum case, non-trivial operator mixing also occurs beyond leading order in lattice perturbation theory. The operator-mixing of the three different continuum operators . . . c into a single lattice operator . . . L follows from Beyond the operator mixing that also involves the i EB -contribution, the convergence of the BB -operator to the continuum was accelerated compared to [16] by replacing the clover-clover BB -correlator with a single-plaquette expectation value. Although this decreases statistical power (4 instead of 16 combinations of single plaquettes), it renders all four measured correlators at the correct lattice spacing, there are no more contaminations of operators at the wrong separation. With this in hand, we were able to also achieve a valid continuum limit for BB at g 2 3d L = 0.25. The discussion of different sources for errors in [16] also applies here.

Lattice vs. perturbation theory
At T = 100 GeV and N f = 5, we expect to be deeply in the perturbative regime. Therefore, this serves as a good starting point for comparing our perturbative results to the updated lattice data in tab. 1. Figure 2 (left) shows that our EE lattice measurements agree with the NLO analytical result within few multiples of the error up to separations of g 2 3d L = 2.0. Although the LO estimate was already quite close to the analytical data, the NLO shifted it by a small margin on top of the non-perturbative solution. The relative correction induced by the NLO is small, so it seems reasonable to assume good convergence of the perturbative series at small g 2 3d L. The same holds for BB in fig. 2 (right). In fact, the agreement of the lattice data with the perturbative result is surprisingly good, as one would expect contributions from the generically non-perturbative magnetic sector to affect BB the most of the three correlators. 9 To avoid confusion, we updated the values of tab. 2 in the arXiv-version of [16], except from the i EB measurement and the g 2 3d L = 0.25 data point of BB . Note also that this publication uses slightly different conventions for the normalization of the fields, such that Z 3d g has now mass dimension 1 instead of 2 in [16].   Table 1. Continuum-extrapolated results for the correlators EE /g 6 3d , BB /g 6 3d , and i EB /g 6 3d , at four tuples (T, N f ) of temperature and number of massless quark flavors over a range of separations g 2 3d L. (Z 3d g ) latt sums the integrated correlators according to (3.4) with small-g 2 3d L cutoff for EE and BB at g 2 3d L min = 0.25 and for i EB at g 2 3d L min = 0.5. Note that the scalar field Φ on the lattice is rescaled such that it is dimensionless. Therefore, the three correlators on the lattice -and Z 3d g accordingly -feature slightly different normalizations than in perturbation theory (cf. (3.4), app. B). (Z 3d g ) cutoff contains the entire integration of lattice data points, LO subtraction and long-g 2 3d L-tail integration in the range g 2 3d L > g 2 3d L min , whereas (Z 3d g ) full ranges over the full integration domain (cf. sec. 5.3). The subtraction of subleading O(dL L/L 2 ) divergences still contains a logarithmic dependence on the IR cutoff at the respective g 2 3d L min , which will be removed when supplied with four-dimensional results. The data for Z 3d g have three errors: the Monte-Carlo integration causes an error for the integration range in which lattice points are available (first bracket), the error for the integration of the tail caused by the Monte-Carlo error of the points to which the functional form is fitted (second bracket), and finally the error introduced by the finite difference integration via the trapezoidal rule (third bracket). The correlator i EB vanishes at leading order. Therefore an analysis of its convergence would require a NNLO result, which is not available as of now. Figure 3 collects lattice data at all four pairs of (T, N f ) and includes their corresponding NLO predictions. It can be seen that with smaller temperatures -and larger coupling -the onset of perturbative behavior in the UV decreases to smaller g 2 3d L. Figure 3 also shows that perturbation theory qualitatively predicts well i EB at small g 2 3d L and high temperatures where the coupling is small due to asymptotic freedom. With lower temperatures and larger separations the agreement with our lattice data becomes gradually worse until perturbation theory even fails to predict the correct sign of i EB at small separations and the smallest two temperatures. What saves the day is that on the one hand, i EB is numerically suppressed compared to EE and BB , so its overall impact on Z 3d g is not large, as seen from fig. 3(lower right). On the other hand, the tree-level contribution for i EB vanishes. Therefore, we evaluate only one non-trivial order for i EB and cannot analyze how well the perturbative series converges or assess the quality of our perturbative estimate for i EB . Bearing this reasoning in mind, we accept the mismatch of the lattice data for i EB and its perturbative prediction in fig. 3.
To mitigate the discretization effects of the numerical integration via the trapezoid rule, we add and subtract our NLO result to the integration dLL (corr) latt = dLL (corr) latt − (corr) NLO + dLL (corr) NLO , (5.2) with the NLO form of the correlators in (4.12), (4.13), and (4.14). The last integral can be done analytically, whereas the difference in the first integral is numerically much better behaved: O g 4 3d /L instead of O 1/L 3 for EE and BB , or g 2 3d /L 2 for i EB , respectively. At this point, we would like to caution the reader that this procedure is purely motivated to accelerate numerical convergence and has no physical meaning, in contrast to the subtractions in the following subsection.
As elaborated in [16], it is necessary to model the large-g 2 3d L tail of the correlators in order to perform the dL L integration up to ∞. For EE and BB , the functional form, motivated by [30], is A with the fitting constants A and B. Considering i EB , we find that the data rather follows with the respective fitting constants A ′ and B ′ . As already argued above, the impact of i EB on Z g is small. Also, the error associated with the large-g 2 3d L tail given in tab. 1 Table 2. Definition of (Z 3d g ) UV , as given in eq. (5.5). The second and third columns give the integrands for the three contributions -which are to be added as per eq. (3.4) -in different intervals of g 2 3d L and with IR cutoff g 2 3d L min . We recall that at leading order BB LO agrees with its subtraction and i EB LO vanishes. Moreover, BB subtr NLO and i EB subtr NLO vanish since the respective correlators are UV-finite.
amounts to roughly 100% of the contribution of the large-g 2 3d L tail. Therefore, one can exercise a bit of freedom in the choice of the large-g 2 3d L tail modeling. Moreover, we found that switching between an effective asymptotic description and our lattice data was best performed at g 2 3d L = 3.0 for EE and BB , but for the fit of i EB it was most stable when neglecting the noisy g 2 3d L = 3.0 points and switching already at g 2 3d L = 2.5.

Subtraction scheme
As explained in sec. 4, we need to subtract the UV-divergent limiting behavior from the EQCD calculation. That is because, in dealing with the contribution of the scale T at O(g 0 ) and O(g 2 ) (which we have not computed yet), it is impractical to separate the contribution of the Matsubara zero mode from that of the other modes. Hence, that T -scale contribution ends up containing a zero-mode contribution computed without any resummation of the screening mass: it thus must agree with the m D → 0 limit of the EQCD calculation. Hence, this m D → 0 or UV limit of the resummed EQCD calculation must be subtracted from our non-perturbative evaluation, so that the latter could, in our upcoming paper, be summed with the contribution of the scale T at order O(g 2 ). This subtraction poses no problem for the leading 1/L 3 UV behavior in eq. (4.6), which can be safely integrated in dL L up to infinity. As we argued, that is not the case for the subleading 1/L 2 terms in eq. (4.16), which are IR divergent. We introduce an IR cutoff for the NLO subtraction, to be absorbed once the T -scale contribution becomes available. For convenience, we set this cutoff to the respective smallest lattice separation g 2 3d L min of each correlator, i.e. 0.25 for EE and BB , and 0.5 for i EB ; see tab. 1. Therefore, we modify the values for the correlators in tab. 1 by subtracting the leading UV counterterm (4.6) over the entire range, and add the NLO perturbative result, subtracted by the NLO divergence, at g 2 3d L < g 2 3d L min , viz.
The precise form of the UV completions (in EQCD, i.e. unphysical) (Z 3d g ) UV is given in tab. 2. Bear in mind that the LO subtractions EE subtr LO and BB subtr LO are already factored into Z 3d g cutoff up to the UV cutoff g 2 3d L min . The integrated UV completions can be found in tab. 3. The final results of this paper, the gauge condensate with all necessary subtractions of unresummed zero-mode-contributions, and the IR divergence induced by the 1/L 2 -subtraction regulated by the respective g 2 3d L min , are displayed as (Z 3d g ) full in tab. 1.
Revisiting tab. 1 after these subtractions, we discover that all values of Z 3d g full except the T = 250 MeV one are consistent with zero within errors. This stresses once more the good quality of our perturbative result and the small deviations of our lattice result from the perturbative curves in the part of the integration range that dominates Z 3d g . The final reasoning and intuition about Z 3d g will be postponed to a forthcoming publication where the four-dimensional perturbative calculation is carried out [31].

Discussion
The correlator of two insertions of the color-force is a crucial ingredient for the computation of asymptotic jet masses. In particular, it receives non-perturbative contributions from the magnetic and electrostatic sector of hot thermal QCD. In the present paper, we provided a twofold calculation of this correlator in EQCD: an analytical calculation to next-to-leading order as well as a non-perturbative simulation of all involved correlators.
The NLO calculation is needed for the final goal of integrating the force-force correlator from zero to infinite separation. On the one hand, it matches well at small separations with the non-perturbative lattice data for the numerically dominant contributions. Therefore, the NLO result improves the modelling of the small g 2 3d L region, which is inaccessible to lattice simulations. On the other hand, it allows to understand the nature of UV behavior of the EQCD result. On dimensional grounds, the NLO in EQCD is the last order which can give rise to UV divergent terms in Z 3d g . These divergences need to be subtracted to avoid doublecounting of degrees of freedom in EQCD and in the outstanding parts of full QCD. Ultimately, our goal is to subtract the UV limit from the analytical EQCD calculation, since we know that EQCD is unphysical in the UV, and replace it with the correct UV behavior obtained from full (perturbative) QCD.
We found that presenting the complete EQCD result seemed a natural part of the required work, which has to be finally supplied with the aforementioned perturbative calculation of the contribution of the scale T at O(g 2 ) in full QCD to deliver a meaningful non-perturbative result for Z g and consequently m 2 ∞ . In particular, having the full QCD result at hand allows one to remove the unphysical UV behavior of EQCD, which has been recast by the subtraction procedure of sec. 5.3 as a dependence on an intermediate regulator g 2 3d L min .
In an upcoming paper, we will deal with that missing contribution by determining Z g at NLO in a naive perturbative expansion in full QCD, which is appropriate for capturing the contribution of the scale T . We anticipate that practicality concerns will most likely force us to tackle such a computation in momentum space, forgoing the possibility of keeping the L dependence in intermediate steps. Hence, the expected IR divergence would show up at small momenta, and we would need to translate that in the L-space form of sec. 5.3. That would lift the logarithmic, unphysical dependence on g 2 3d L min introduced by eq. (5.5), leading to a finite, scheme-independent result for the asymptotic mass which incorporates the IR EQCD contribution to all orders.
Here K ν (z) is the modified Bessel functions of second kind and coordinates are chosen so that r points in the z-direction with r = |r|.

A.2 Diagram (b)
The simplest NLO contribution is diagram (b) which is depicted and calculated as follows: using the free propagator integrals (A.1) and (A.2) and the exponential integral function with respective limits: with E 1 (x) = −Ei(−x) for x > 0. In both contributions the corresponding LO operator interaction of diagram (a) factorizes from the internal propagator between the two Wilson line insertions. As exemplified in the internal z 1,2 integration of (b) EE in eq. (A.3), a UV divergence of individual integrals is of opposite sign which renders the overall EE and BB contributions separately finite.

A.3 Diagram (c)
As depicted in fig. 1, diagram (c) requires the self-energies of the adjoint scalar Φ a and the spatial gauge bosons A a i ; diagrammatically given by where a curly line is a spatial gauge boson, a solid line an adjoint scalar, and a dotted directed line a ghost field. The analytical forms of the self-energies are where the gauge boson self-energy, Π A a , is transverse in three dimensions. Addtionally, we set d = 3 − 2ǫ and employ the master integral [29] The gauge boson self-energy agrees with [19] in the limit of soft external momenta p 2 The remaining integrals needed for the p-integration of eqs. (A.9) and (A.11) follow by inserting eq. (A.12) in the following , (A.17) (1 + mr) ln(2mr) + γ E + (1 − mr)e 2mr E 1 (2mr) , (A.20) which we evaluated using the residue theorem along the integration contour as in [32]. For the last integral I 2;0mm only its finite second r-derivative is needed. As a result, we obtain for the diagonal contributions in Feynman gauge: To understand their UV asymptotics we expand for small m D L → 0 which vanishes for (c) BB in (A.26).

A.5 Diagram (e)
Conveniently diagram (e) is finite. The scalar contribution originates from the non-Abelian part of D x Φ and the gauge contribution originates from the non-Abelian part of the fieldstrength-tensor F a xz ⊃ g 3d f abc A b x A c z . Their momentum space integral representations are where we implicitly used that the gluon propagator is spatially diagonal in Feynman gauge. The result is logarithmically divergent which, however, vanishes after considering opposing signs of EE and BB in (3.4).

A.6 Diagrams (f ), (h)
Diagram (f ) corresponds to the Fourier transform of the EQCD three-point vertices. Starting from the electric scalar contribution in momentum space integral representation, the strategy is to first integrate over the z-coordinate of the single A z field sourced from the Wilson line (A. 30) In obtaining the final line, we have exploited the p ↔ k symmetry of the integrand to reshuffle the e ikzL term into an e ipzL . As the integrand was, prior to this step, well behaved for p z = k z , we can treat the apparent pole at p z − k z introduced by this reshuffling with a principal value prescription P (cf. [33]). To perform the momentum integrations in eq. (A.30), we proceed as follows. We (i) combine the two k-dependent denominators with a Feynman parameter x, (ii) perform the d 2 k ⊥ and the P-regulated dk z integrals, followed by the one over the Feynman parameter x which yields We obtain Diagram (h) EE and its mirrored partner vanish because each transverse-momentum integration has an odd integrand, The magnetic contribution (f ) BB is strategically analogous to its electric counterpart. The z-integration over the single A z drawn from the Wilson line yields The corresponding master integrals of the Fourier transform are The sum over the corresponding terms in eq. (A.34) is not finite, i.e.
2 × (f ) BB = −g 2 3d C 2 A C F I 00 011 + I 00 101 + I 00 110 + 2(d − 2)I 11 111 − 2 I 0;000 where the final term on the first line is the contribution of the simpler, final two terms on eq. (A.34) and its master is the massless I 0;000 = 1/(4πL) 2 from eq. (A.17). As we now show, this UV divergence is cancelled by an opposite one from diagram (h) BB , which reads so that the combined result is finite again At NLO two contributions to i EB of mixed electric-magnetic correlations arise. One of them is (f ) EB and its mirrored partner (f ) BE . Their result in momentum space is integrated over the z ′ -coordinate of the Φ field sourced by the Wilson line (A.43) The first three terms can be evaluated by using a bare (m D → 0) subtraction scheme where the first bracketed term is evaluated in d = 3 and the second in d = 3 − 2ǫ dimensional regularization. The last term in eq. (A.43) splits into a part treatable with the same scheme and an additional finite part proportional to with the hyperbolic sine integral function Shi(x) defined below eq. (4.14). To evaluate this finite contribution, we went to position space and inspected its asymptotic behavior for small and large values of m D L. Thus, we could extract the analytic form of eq. (A.45) which agrees with the numerically integrated result. After summing all terms in eq. (A.43), we obtain The second mixed i EB contribution is (h) EB which has the momentum space representation after the z-integration of the single Wilson-line sourced field Its scale dependence and ǫ-poles are compensated by diagram (f ) EB (A.46). We collect here the sum of all EB-terms, which is where [BE] collects mirrored contributions.

A.7 Diagram (g)
The scalar electric contribution (g) EE gives rise to the following Fourier transform including the evaluation of the master integrals The magnetic contribution (g) BB vanishes in contrast to its (g) EE -counterpart. This is seen explicitly in where the last line holds in dimensional regularization and vanishes for d = 3 since the integrand is finite. It was derived using integration-by-parts identities at two-loop level similar to eq. (A.15). The relevant relation is (A.52)

B Simulation details
In order to carry out non-perturbative computations, we discretize the continuum action of EQCD (3.1) on a three-dimensional numerical grid with lattice spacing a, S EQCD,L = β x,i>j 1 − 1 3 x,ij + 2 x,i where U i (x) is the gauge link in spatial direction i, connecting lattice sites x and x + aî. We rescaled the adjoint scalar field to its lattice version Φ L such that its wave-function normalization Z Φ = 1 is always enforced. The subscript L will be dropped in the context of lattice calculations. Analytical calculations in EQCD lattice perturbation theory yield expressions for the counterterm of the inverse gauge coupling β, the quartic coupling δx, and the multiplicative mass (Z 2 ) and quartic (Z 4 ) renormalization that compensate for discretization errors up to O(a) [34]. A semi-analytical computation of δy completes the O(a)-improvement at the Lagrangian level [35]. At the operator level, we use the lattice implementation of the modified Wilson line (3.3) described in [36], with the O(a)-improvement delivered in [37], to connect the respective pairs of operators. The 'color-electric' field operator is discretized using gauge-invariant central derivative (see [16]), whereas the 'color-magnetic' field operator is calculated using the clover discretization [38,39]. For the BB -correlator, only the single-plaquette correlations of the pair of clover operators are used, as elaborated in sec. 5. Consequently, the only remaining sources of errors at O(a) are of the form (5.1) and can be determined in an overall fit, as outlined in sec. 5. Our numerical implementation is based on openQCD-1.6 by Martin Lüscher [40]. Using a combined update of one heatbath sweep succeeded by four over-relaxation sweeps through the volume, we update the lattice sites in a checkerboard ordering. We use the multi-level algorithm proposed by Lüscher and Weisz [41] to reduce the noise in relation to the signal for non-local operators in lattice gauge theories. In particular, we divide our volume in four sub-volumes along the z-axis and freeze the surfaces between them for 80 combined heatbath/over-relaxation sweeps, before we allow a single combined sweep through the entire volume. For further details, see [16,35,42]. We give the raw data at finite lattice spacings a in an online repository [43].
We provide the parameters of our simulations in tab. 4. As argued in [44], EQCD possesses a mass gap and therefore finite volume effects are exponentially suppressed. Thus, one can keep finite volume effects under control by maintaining a sufficiently large volume along the rules of thumb given in [44]. We measured all three correlators on the same lattice configurations introducing correlations among different g 2 3d L and correlators that had to be accounted for in the jackknife error analysis by keeping the binned jackknife data until the numerical integration and calculating the error thereafter. Equal computational resources were spent on all four scenarios (T, N f ); the difference in statistics can be explained by a slight increase in the acceptance rate for smaller x, since the quartic self-coupling of the scalars is taken into account via a Metropolis step in the scalar heatbath-update.  Table 4. Parameters for all EQCD multi-level simulations. We give the continuum expressions for the rescaled quartic scalar self-coupling x cont and the rescaled screening mass y cont . A conversion into lattice units is done using the respective lattice-continuum relations [34,35].