The nonperturbative contribution to asymptotic masses

In gauge theories, charged particles obey modified dispersion due to medium interactions (forward scattering), leading at high energies $E \geq T$ to an asymptotic mass-squared $m^2_\infty$. We calculate the infrared part of this mass nonperturbatively for the theory of the strong interactions, QCD, through a lattice treatment of its low-energy effective description, Electrostatic QCD (EQCD). Incorporation of these results into a nonperturbative determination of the effective thermal mass will require a still-incomplete next-to-leading order perturbative matching of this quantity to full QCD.


Introduction
The Quark-Gluon-Plasma (QGP) is an exotic state of matter created at collider facilities [1][2][3][4] and delivering important insights into the nature of strongly interacting matter. The experimental effort to measure the features of the QGP with the highest possible precision goes hand in hand with the theoretical goal of giving precise predictions for these very measurements. The properties of the QGP can in essence be divided into two classes: thermodynamics and transport properties. Typical examples for the former are the pressure, perturbatively determined up to O(g 6 ln g) [5][6][7][8][9][10] and investigated nonperturbatively on the lattice [11,12], and charge susceptibilities and other quark-number correlations, which are well investigated on the lattice [13]. Since thermal equilibrium is time-translation invariant, these results can be derived in Euclidean time, which simplifies the calculations considerably. In contrast, the calculation of transport properties usually requires an explicit treatment of Minkowski time, which is conceptually more demanding. Nevertheless, there has been significant progress in recent years, for instance for the shear viscosity η s [14,15] or the thermal photon production rate [16,17].
Anther interesting transport phenomenon to investigate is how the medium modifies a jet, whose theoretical treatment is for example summarized in [18]. For our purpose, jets are strongly interacting particles that carry a large momentum and move at (almost) the speed of light. Because of the large momentum scale, asymptotic freedom of QCD suggests that the strong coupling constant g is small [19,20], and that they can be well described by an expansion in powers of g. In traversing through the medium, however, they receive highly nontrivial modifications [21,22], pictured by a large number of soft bumps between jet and medium constituents. This interaction eventually leads to jet broadening, i.e. a broadening of the jet's transversal momentum distribution. Related to that, jet quenching refers to a jet losing energy to the medium [23,24], predominantly via stimulated emission due to multiple soft scatterings with medium constituents.
A quantity that influences both phenomena is the asymptotic mass. It was already found by Klimov and by Weldon that massless particles, be it gluons [25,26] or fermions [27,28], tend to follow the dispersion relation of massive particles with their mass m 2 p depending on the momentum p and generated by the interaction of the particle with the medium (essentially, forward scattering). At large momenta, the particle's velocity in the plasma frame approaches the speed of light v → (1,p), and this thermal mass approaches a constant which we will call the asymptotic mass m 2 ∞ . At the leading perturbative order, m 2 ∞ can be computed from the imaginary part of the transverse Hard-Thermal-Loop-resummed self-energy Π T near the light cone [29,30]. The gluon's asymptotic mass, for instance, reads [31]  to first perturbative approximation. Here m 2 D = (2N c + N f )g 2 T 2 /6 is the Debye screening mass squared.
The expression (1.2) can be generalized to a particle in an arbitrary representation R and related to a combination of two condensates [32,33], a gauge contribution Z g and a fermionic contribution Z f : with the (running) gauge coupling g, and the Casimir operator of the representation C R of the particle that makes up the jet. The two condensates are [32,33] c − 1 is the dimension of the adjoint representation, in which the gauge bosons live, and d R is the dimension of the representation R of the fermion, for instance the fundamental representation of SU (3) for quarks in QCD, and the angular brackets represent an average over all thermal fluctuations. Note that our definition of the QCD field strength F σµ is the geometrical one, F µν = i[D µ , D ν ] with D µ = ∂ µ − iA µ ; g 2 is absorbed into the coefficient on the field strength in the action, S = 1 2g 2 TrF 2 +. . .. These condensates represent the contribution of forward scattering in the medium, with 1/(v · D) or (v · D) 2 accounting for the lightlike (eikonalized) propagation through the medium and ψ, ψ or v σ F σµ representing an interaction with a medium excitation. The leading-order, long-distance approximation to these condensates reproduce the forward-scattering part of the hard thermal loop effective action along the light cone.
The gauge part of (1.3) can be related to an integral over the correlator of two covariant Lorentz force vectors in position space [15] where the x + integration and the modified adjoint Wilson line U ab A arise from 1/(v · D) 2 . Locations in space-time are described by three-vectors (x + , x − , x ⊥ ) in light-cone coordinates. Because the jet is highly relativistic, we have x − = t − z = 0 in the plasma frame. The problem therefore effectively reduces to three dimensions, two in the transverse plane x ⊥ , where rotational invariance is present, and a coordinate that encompasses the elapsed time t for the jet as well as its covered distance z in the medium x + = (t + z)/2 ≡ L. We use the convention that the z-component of k is the jet's direction of propagation, so v = (1, 0, 0, 1). At this point, we would like to briefly build up some intuition for the specific form of (1.5). We can distinguish two sorts of jet-medium interactions. First, there is true scattering, where the transverse momentum of the jet particle changes; this is described by the transverse momentum broadening kernel C(b ⊥ ) [22]. But there is also forward scattering, where the jet particle temporarily changes direction before scattering back into its original state (or more technically, a scattering creates an amplitude in a different momentum state, which subsequently scatters again to return to the original momentum, introducing a phase shift). Such forward scattering is what causes a dispersion correction. It can be understood qualitatively as the jet particle making a scattering-induced detour, as pictured in Fig. 1. The sum of all such forward scattering possibilities induces an effective mass term [23].
Relativistic kinematics tells us the length of the detour where m is a dummy jet mass that will eventually drop out of the calculation, γ(v) is the gamma-factor of the jet at velocity v, and T is the time it takes for the jet to cover the distance L. If we consider this problem from a frame of reference in which we do not take into account the additional detour, the kinematic invariant reads 1 in the ultra-relativistic limit, in which also the contribution of longitudinal forces to the effective mass vanish. Consequently, the (transversal) gauge force must appear twice in the correlator resulting in the mass. Furthermore, a full quantum treatment requires a sum over all possible paths instead of just considering one as in Fig. 1. Hence, our estimate for the effective mass yields reflecting the qualitative features of (1.5).
The state-of-the-art perturbative result for Z g at next-to-leading-order (NLO) [33] is The Debye mass [26], in turn, depends on the coupling g: in QCD (N c = 3) it reads (1.10) Formally the NLO correction is therefore suppressed by a single power of the coupling g. Numerically, it is not suppressed even at T = 100 GeV temperatures, indicating a very poor convergence for the perturbative series. This is a familiar situation for perturbative calculations in hot QCD, and as usual, it arises because there are large corrections arising from the large-occupancy, infrared sector of QCD. Specifically, the convergence of the perturbative series is spoiled by the high occupancy of the gluon 0-Matsubara-frequency mode. When this mode circulates in a loop, it introduces an additional factor of n B ∼ 1/g [34][35][36]. A possible solution to this problem is provided by treating the badly-convergent contributions separately in a framework called Electrostatic Quantum Chromodynamics (EQCD) [37]. In contrast, the situation for Z f appears to be better under control. Fermions do not feature a 0-mode so the above reasoning does not apply here. Indeed, it was found that [33] so the assumption of a convergent series for the fermionic part of m 2 ∞ may be justified. According to [38], correlators like (1.4) can be computed in a setup that quantizes the plasma on the light-cone. Their dynamics is dominated by multiple soft scatterings that have to be resummed. But because they are infrared, they can be computed in three dimensional EQCD instead of four dimensional full QCD in real-time, with modified Wilson lines that account for the parallel transport of the gauge configurations. In the past, this technique has been predominantly used for the computation of the transverse collision kernel C(b ⊥ ) [38][39][40]. Similarly to the calculation of C(b ⊥ ), the EQCD operator agrees with its full-QCD-equivalent only in the infrared (IR) regime where EQCD is supposed to be a valid description of full QCD. In the ultraviolet regime, however, a difference between the two operators will emerge. This difference should be under perturbative control and can be included as part of a matching calculation between thermal QCD and EQCD, when the matching is extended to include the operators of Eq. (1.5).
In the course of the present work, we will evaluate the soft operator in EQCD and leave the purely perturbative matching calculation for later investigation. As for C(b ⊥ ), perturbation theory, even in EQCD, is not always sufficient to deal with the nonperturbative nature of some high-temperature phenomena. Sometimes, it is necessary to solve EQCD directly on the lattice, which has been done recently for C(b ⊥ ) [39,40]. This will also be our method of choice for the determination of the soft part of asymptotic jet masses in the following.
We address the problem as follows. In Sec. 2, we outline how we would compute Z g if there were no interactions. To this end, we briefly introduce how our setup translates into EQCD, show how the observables appear in this theory and finally present the free solutions. Solving the full interacting theory requires discretizing the theory and the corresponding observables on the lattice, done in Section 3. This part of the present work also outlines how to deal with the contributions of Eq. (1.5) at the shortest and longest distances, where the lattice determination becomes problematic. Section 4 involves a detailed presentation of our results together with a thorough discussion and, where possible, quantification of the errors that our results are fraught with. The consequences of our results and how to improve on them is discussed in Sec. 5. Appendix A provides a detailed overview of our lattice configurations.

Electrostatic QCD
Electrostatic Quantum Chromodynamics, or EQCD for short, is an effective field theory for hot QCD. At high temperatures, the scaling behavior changes, for lengths longer than 1/T or energies lower than T , to the behavior of a 3D theory, which is more strongly coupled in the IR. Within the Matsubara formalism for describing thermal QCD, this arises because the usual frequency integral dω/2π is replaced by a discrete sum T ω=2nπT , with the n = 0 (zero Matsubara frequency) mode massless (at tree level), leading to an O(2πT /p) enhancement of small momentum p parts of diagrams. For the electric sector this is cut off by the Debye mass, but for the magnetic sector it is not cut off until magnetic degrees of freedom become strongly coupled and confine. For an overview see Refs [34][35][36]. The idea of EQCD is to integrate out -in a Wilsonian renormalization-group-sense -all Matsubara modes in QCD except for these 0-modes, which can then be treated nonperturbatively.
Since there is only one Matsubara-mode left, one has effectively eliminated the Euclidean time direction. Therefore the transition to EQCD is often called dimensional reduction. A time direction that only consists of one time layer does not allow for derivatives in time direction any more. Consequently, this direction no longer features gauge symmetry, which is what allows the temporal gauge field A 0 , now responsible for all electrical phenomena, to acquire a mass m 2 D . Loop level effects induce a self-coupling for this field, denoted as λ. A scale is set by the now dimensionful gauge coupling g 2 3d = g 2 T . The full action of EQCD in the continuum reads The parameters in (2.1) can be related to the full theory via a perturbative matching procedure [37]. We employ a O(g 4 )-accurate result [41]. Input parameters are the temperature T and the number of dynamical massless fermion flavors N f translating into the pair of (dimensionless) EQCD parameters x ≡ λ/g 2 3d and y ≡ m 2 . Note that m 2 D varies logarithmically with scale due to 2-loop effects; therefore in defining y we must specify the renormalization point. Following convention we choose µ = g 2 3d . We will consider the values given in Tab. 1. EQCD also has a nontrivial phase structure [42,43], with a transition between a Z 3 -symmetric and broken phase that is of second order at large x and turns into a first order one for smaller x, including all physically relevant x values. Note the (x, y)-pairs which describe full QCD prove to lie in the region where the 3D theory prefers the Z 3 -broken phase. Therefore we are to consider metastable Z 3 -symmetric states in a parameter range where the broken phase is actually thermodynamically preferred. Thus choosing sufficiently large volumes is crucial; they keep the system from accidentally falling into the broken phase. It is believed that dimensional reduction to EQCD is valid for all T 2T c and possibly even further down [44], where T c is the transition temperature of full QCD.  Table 1. EQCD parameters at four different temperatures and number of fermion flavors.

Observables
As already mentioned, one promotes (1.5) to EQCD by replacing Wilson lines along the light cone with their modified EQCD-counterparts [38] and F i0 with iD i Φ 2 . Applying rotational invariance in the transversal plane, we find Using the modified Wilson loop Eq. (2.2) comes with the advantage of the correlator being free of logarithmically UV divergent perimeter law suppression which is cancelled by the scalar insertions Φ a T a in Eq. (2.2). However, this correlator by itself is not free of L → 0 UV divergences which have to be cancelled by a suitable subtraction. Since the correlators on the lattice can be measured at only finite L, anyway, we do not perform such a subtraction in this work and leave it to future investigations. As a matter of choice, we work in the fundamental representation of the correlators, related to the adjoint ones via and for the cross-correlators, respectively. In the following, we will drop the tilde on top of the modified Wilson linesŨ for the sake of readability, and write half the right-hand side of (2.4) as Tr E(L)U E(0)U −1 and half the right-hand side of (2.5) as Tr B(L)U B(0)U −1 .
Let us stress again that the argument of Eq. (2.3) agrees with that of Eq. (1.5) only up to IR-safe terms which can and must be determined in a matching calculation between full thermal QCD and EQCD. This calculation has not been carried out for these specific operators, though it should proceed along similar lines to the matching already conducted for C(b ⊥ ) in [15,38].
Through insertions of the scale g 2 3d , everything can be recombined into dimensionless quantities and we end up with the master formula for our observable but we will omit the commutator for the sake of readability 3 In an earlier version of this publication, we were not aware of the cross-correlations between DxΦ and F xz . They do not appear at tree-level; due to the Φ → −Φ symmetry in EQCD, an odd number of Φ-field has to be sourced by the Wilson line. We thank Jacopo Ghiglieri for pointing that out to us.
where g 2 3d L EE min , g 2 3d L BB min , and g 2 3d L EB min are minimal separations at which we can measure the respective observable in the EQCD effective picture on the lattice. Beyond that, we have to rely on analytical solutions, as we will elaborate in Sec. 3. We will only calculate the Tr E(L)U E(0)U −1 and Tr B(L)U B(0)U −1 correlators in this work. These numerically dominate over the cross-correlations, anyway.

Free solutions
The short distance behavior of the full, continuum-extrapolated correlators in (2.4) and (2.5) is expected to match the perturbative expectations and eventually, at sufficiently small separations, the free solutions. Since they not only provide an important crosscheck for our data, but also will be of practical use later, we will briefly present their instructive but simple derivation in the following.
We assume the transverse direction to be oriented along the x-axis for simplicity. At tree-level, i.e. at O(g 0 3d ), (2.4) reduces to specializing to SU (3) in the last step, and in a similar fashion for (2.5), Plugging these results into (1.5), we find a leading UV divergence of O(L −3 ) which has to be subtracted. A contribution of O(L −2 ) cancels at tree-level, whereas a O(1/L)-contribution persists, which when integrated over L dL leads to a finite short-distance behavior in Z g . The superrenormalizable nature of EQCD means that, at the next loop order, the small-L behavior can be at worst O(L −2 ). However, if there is no cancellation between the electric and magnetic contributions at this order, there could be complications in the short-distance behavior, leading potentially to logarithmic short-distance contributions in Eq. (2.3). If this is the case, their role must be clarified by the matching to the 4D theory.

Lattice implementation
Our lattice implementation agrees with the one in [40], based on a modification of Martin Lüscher's openQCD-1. 6 [45]. Due to [46] and [43], we are able to discretize the EQCD action on a three-dimensional spatial lattice without any errors at O(a). As an update algorithm, we used a composite of 2 heatbath sweeps followed by 4 overrelaxation sweeps through the volume in a checkerboard-order. To improve the signal-to-noise ratio, we applied the multi-level algorithm [47,48], dividing the volume into four subvolumes along the z-axis and updating each subvolume 80 times before a sweep through the entire volume took place.

Observables on the lattice and their continuum extrapolation
For the F xz -insertions in (2.5), we use the standard clover discretization of the field-strength tensor [49,50], with four 1 × 1 clover leaves. The scalar field derivatives D x Φ in (2.4) are discretized as keeping in mind that the rescaling of the Φ-fields to their lattice equivalents yields another factor of 1/a. Note that the clover magnetic field operator is spatially larger than the electric field operator we use; therefore its O(a 2 ) corrections are larger, and we need more lattice spacings of separation between the B operators before continuum behavior sets in than for the E operators. Consequently we cannot pursue as short distances in the Tr B(L)U B(0)U −1 correlators as in the Tr E(L)U E(0)U −1 correlators. All of our observables are tree-level O(a 2 ) accurate. Nevertheless, there are still O(a) discretization effects present generated by loop effects in lattice perturbation theory. These generate multiplicative O(a)-corrections, which we will denote by Z E and Z B in the following. For details about the discretization of the modified Wilson line, we again refer to [40] and [39]. The O(a)-renormalization of an isolated modified Wilson line, or two wellseparated modified Wilson lines, was calculated analytically in [51]. We use the resulting prescription, which is to scale A 0 in Eq.  An analytical calculation of Z E , Z B , and Z P is probably possible but appears to be technically somewhat complex, since the Wilson line operator is an extended object. However, by analyzing what diagrams could give rise to such linear in a effects, we find that the rescalings Z B and Z E do not depend on the length L (provided it is long in lattice units), and the perimeter contribution Z P is strictly linear in L. Furthermore, the EQCD parameters x and y would first enter at a higher loop order than the NLO effect giving rise to O(a) corrections. Therefore the Z E , Z B , and Z P factors are common to all (x, y) values and all separations. They represent only 3 total fitting parameters over all of our data at multiple (x, y) and L values. Such global parameters make the fit more complicated, introducing correlations between the fits to multiple lattice spacings, separations, and (x, y) pairs. This made the fitting procedure somewhat unstable, which we cured by using a finite-difference solver to minimize χ 2 . The correlations between all datapoints must also be considered when performing the final numerical integrations. But because we only lose fitting 3 degrees of freedom over all data, the impact on the overall error budget is relatively modest.

Integration of lattice data
Even without the contributions from the 4D to 3D matching calculation, we can still attempt to perform the integral shown in Eq. (2.6), at least for L above some small-L cutoff. Our lattice data is available only within a finite range of separations g 2 3d L. We integrate this data using the trapezoidal rule. Ideally we would want data at a densely spaced set of distances L, but this is impractical because the L spacing is limited by the lattice spacings of the coarsest lattices used. Also, using denser-spaced L values would not actually help very much, because the different L-value data are from the same simulations, and data at nearby separations L are highly correlated. Nevertheless, the finite spacing in g 2 3d L introduces another error in Z g which we have to bear in mind. We A. We will eventually find that using the trapezoid rule only introduces a subdominant discrete-integration-error.

Integration of the tails
The correlator (1.5) requires an integration of the condensate up to infinite separations. Since the signal-to-noise ratio drops dramatically with increasing L, correlators at separations beyond g 2 3d L ≥ 3.0 are not determined robustly. Nevertheless, knowing the analytical form of the correlators at large separations L, we can fit the tail to the largest data we have available. Similar things have been done, for instance, in the context of the hadronic contribution to the muon g − 2-factor [52]. However, their considerations cannot be precisely applied in our case since we do not know the exact form of the higher excited states. For the magnetic contribution, Laine and Philipsen [53] conjecture a functional form Even though they consider a conventional Wilson line connecting the two B-fields, this form is still applicable. However in our case the coefficient B has a valid continuum limit, rather than having logarithmically UV divergent perimeter contributions. Each coefficient is nonperturbative; we cannot predict them. Based on similar physical reasoning, and consistent with the free solution (2.7), we expect the functional form of the large-separation limit of Tr E(L)U E(0)U −1 also to be (3.4), with different coefficients to be determined from the fit. We choose our fitting range for the large-g 2 3d L-tail to be [1.5, 3.0]. This is hard to justify a priori, since the correlation length above which the lower bound of the interval must be is precisely ξ ∼ 1/m ∞ . However, we will eventually find that 1.5 g 2 3d /m ∞ . At the other end, we do not have data at separations below certain minimal separations g 2 3d L E min and g 2 3d L B min . Below these minimal separations, one would have to connect to an EQCD perturbative expression. Since the (yet unknown) difference between the EQCD effective description and the "true" full-QCD version of (1.5) is expected to contribute significantly there, we leave the consideration of this region entirely for future study.

Results
Before we give our numerical results, we present a thorough discussion of the multiple possible sources of both statistical and systematic errors throughout our calculation. As numerical results, we provide continuum-extrapolated data for Tr E(L)U E(0)U −1 and Tr B(L)U B(0)U −1 in Figs. 3 and 4, the corresponding values can be found in Tab. 2, and details about the simulations appear in Tab. 3 of App. A.

Error estimation
The estimation of the error has to take several sources into account. They are: We will now discuss these aspects in detail, and give and explain quantitative estimates where possible, see Tab. 2.
The Monte-Carlo error of Tr E(L)U E(0)U −1 and Tr B(L)U B(0)U −1 can be determined by a standard binning and Jackknife analysis. Since the continuum extrapolation was performed in a grand fit, one has to account for correlations of all lattices at all separations g 2 3d L. Even if their g 2 3d a naively would not contribute, it does so via its influence on the values of Z P , Z B , and Z E . This error, ∆ MC , is given in the first brackets in Tab. 2.
Another dominantly statistical error is created when we match the tails of our models for Tr E(L)U E(0)U −1 and Tr B(L)U B(0)U −1 , see (3.4). To be precise, the origin of this error is twofold; the assumption of the model itself introduces a systematic error, and the determination of the model parameters from our data at large g 2 3d L generates a statistical error. A valid estimate for this is provided by determining the coefficients in (3.4) by a fit, performing the integral, repeating this procedure for each Jackknife bin, and running the standard Jackknife-analysis over the different values for the tail-integral. This procedure includes the correlations between the fitting parameters in Eq. (3.4). In practice, we found that the statistical error amounts to about 100% of the correction induced by the large-g 2 3d L-tail, which seems to be a reasonable estimate of the overall error of the large-g 2 3d L-integration. This error, ∆ LT , is given in the second brackets of Table 2. Our data in the available window was integrated using the trapezoidal rule. The remainder ∆ TR of such a discretized integration bounded by g 2 3d L min and g 2 3d L max reads where the maximum of the second derivative of the integrand I in (2.3) within the interval g 2 3d L min , g 2 3d L max has to be found, involving the additional factor of g 2 3d L from the L dLintegration. This is tricky since the overall functional form of our data is precisely what we do not know. However, we know the functional form of the tails which we take as a basis. This might not be a rigorous choice but we will see that this error is in any case subdominant compared to other error sources. For this consideration, we extend the range of validity of the free solutions up to g 2 3d L = 1.0 and the range of validity of the large-g 2 3d L-tails down to the same value. Consequently, we have three trapezoids of length δ(g 2 3d L) = 0.25 with free-solution estimates of the second derivative and four trapezoids of length δ(g 2 3d L) = 0.5 with tail estimates. We give this systematic error in the third brackets in Tab. 2.
Furthermore, the dimensional reduction to EQCD also introduces errors. These errors are of two types. First, there is the precision with which the matching is computed for the Lagrangian of EQCD. These corrections are formally suppressed by O(g 4 ) relative to the leading behavior. Second, there is the accuracy of the matching for the specific operator we are considering. So far this has only been performed at tree level, which as we have already noted is insufficient. A calculation at the next order would represent parametrically O(g 4 ) contributions to Z g , which is formally O(g 2 ) relative to the leading behavior. Without this matching there are formally O(g 2 ) relative errors from the matching calculation; with an NLO calculation this would be reduced to formally O(g 4 ) errors. The reason an NLO matching is needed is that nonperturbative EQCD effects are also formally of relative order O(g 2 ); therefore without an NLO matching calculation for the operator, we have not improved the formal order of the uncertainty. It is difficult to estimate the size of these (systematic) matching errors. An NLO calculation would shed considerable light on this point.
Last but not least, the current perturbative determination of Z f neglects terms suppressed by a factor of g 2 . These could be determined in an NLO calculation of this condensate within the full 4D theory.
All in all, the values of Z g in Tab. 2 are of the form best estimate(∆ MC )(∆ LT )(∆ TR ) . Here ∆ MC , ∆ LT , and ∆ TR are respectively the errors from Monte-Carlo statistics including the continuum extrapolation, the long-distance tail, and the use of the trapezoid rule. The overall quantified error of our results is dominated by integration of the longdistance tail. This introduces a relative error of ∼ 1%. However, we assume that analytical uncertainties like the perturbative matching of the action and the operator and the perturbative determination of Z f would not allow for higher precision anyway, so an improvement is not mandatory.

Numerical results
The plots Figs not have data for Tr B(L)U B(0)U −1 at g 2 3d L = 0.25 since the discretization effects are more severe for this operator and we did not compute fine enough lattices to perform a valid continuum extrapolation at this distance. At the other end of the range, we fit the functional form (3.4) to the correlators through the last four points. The fit is dominated by the two innermost of these four points since the signal-to-noise-ratio drops dramatically at increasing g 2 3d L. In turn, this phenomenon leads to some of the points in Fig. 4 at g 2 3d L = 3.0 even being negative and therefore not being displayed in the logarithmically scaled plot. As we can see from their values in Tab. 2, they are still consistent with 0 within a few standard deviations.
We see in both plots that the dependence on the EQCD parameters (x and y from Tab.1) diminishes in going to small g 2 3d L. For Tr B(L)U B(0)U −1 in Fig. 4, this becomes evident by looking at the free solution (2.8) which has neither a dependence on x nor on y, and is jointly approximated by all four scenarios. For Fig. 3 displaying Tr E(L)U E(0)U −1 , the situation is less obvious. Expanding the exponential in  therefore becomes negligible at small distances.
In both cases, we observe a strict hierarchy: the correlators Tr E(L)U E(0)U −1 and −Tr B(L)U B(0)U −1 shrink with increasing x (corresponding to decreasing the temperature). We also find that Tr E(L)U E(0)U −1 and −Tr B(L)U B(0)U −1 become closer to each other as one increases the temperature, which is expected due to the screening mass m D increasing with temperature and reducing scalar correlations over large distances.
Adding Tr E(L)U E(0)U −1 and −Tr B(L)U B(0)U −1 and integrating according to (2.6), as discussed in the previous section, one indeed obtains a strictly positive integrand, yielding a positive Z g .
Since data for the cross-correlations and the NLO matching are still missing, we postpone the calculation of Z g to a later publication.

Conclusion and outlook
In the present work, we have shown how to compute the asymptotic masses m 2 ∞ in QCD at high temperature nonperturbatively. We showed that the real-time continuum expression,  3), to compute the crucial nonperturbative infrared contributions. We then evaluated these operators numerically within EQCD on a 3D lattice, performing a precise extrapolation to the continuum limit. Our data are presented in Table 2, and we discussed a procedure for numerically integrating these data to determine the IR contribution to the asymptotic masses, and for assessing all sources of statistical and systematic error in using the EQCD data.
Unfortunately, we find that to use the numerical results properly, we need to complete the matching calculation for the operators of Eq. (1.5) between full QCD and EQCD a nextto-leading order in the perturbative expansion. Therefore we are in the unusual situation that the numerical study of EQCD is currently ahead of the analytical understanding of the matching calculation. Such a matching calculation would improve the short-distance part of the EQCD calculation by providing an NLO result for the small-separation correlation functions, and it would ensure that all effects up to the desired m 2 ∞ ∼ g 4 T 2 level of accuracy are accounted for. In terms of future improvements, this is clearly the most pressing, since without this matching calculation our results cannot yet be used. Further improvements to our numerical calculation are possible, in particular by sharpening the error bars and computing down to shorter distances, but generally the data appears to be in good shape and the matching calculation is the most urgent need.
Together with the transverse collision kernel C(b ⊥ ) from [43], the then-complete-m 2 ∞ will hopefully contribute significantly to understanding how a jet is modified by the medium that it traverses. We are eager to see the impact of both results on phenomenological transport calculations and how they match experimental observations.  Table 3. Parameters for all EQCD multi-level simulations.