Transverse momentum broadening from the lattice

We present a calculation of the transverse momentum broadening of a highenergy fundamental-representation particle, as calculated nonperturbatively within EQCD using the technique proposed by Caron-Huot and pioneered by Panero, Schäfer, and Rummukainen. Our results are continuum extrapolated and provided at four temperatures: 250 MeV, 500 MeV, 1 GeV, and 100 GeV.


Introduction
Computing transport coefficients of the Quark-Gluon Plasma (QGP) directly from Quantum Chromodynamics (QCD) remains a major challenge in theoretical particle physics. This is mainly due to transport being an intrinsically real-time phenomenon while most of the methods rely on the Euclidean time formalism. Transport coefficients of interest are for instance the shear viscosity over entropy density ratio η s [1] or the jet broadening coefficientq [2]. At high temperatures, one expects the coupling to become small and QCD to become perturbative due to asymptotic freedom [3,4]. However, it has been found that this effect is compensated by highly occupied, and therefore strongly interacting, infrared fields [5].
A major effort was undertaken to compute transport coefficients of QCD plasmas perturbatively at leading order [6,7] and, more recently, at next-to-leading order [8,9]. All computations of next-to-leading order transport phenomena commonly rely on the transverse collision kernel C(q ⊥ ), which represents the rate per unit time for a particle to undergo a scattering which changes the transverse momentum by q ⊥ . The total rate of transverse scattering, and the related momentum-broadening parameter, are (1.1) One can Fourier transform C(q ⊥ ) into C(b ⊥ ), with b ⊥ representing transverse distance (the impact parameter). Casalderry-Solana and Teaney showed that C(b ⊥ ) is determined by the falloff with length L of the log-trace of a Wilson loop with two edges of length ±b ⊥ and two lightlike edges of spatial length L [10]. Remarkably, this clearly Minkowski definition still allows a calculation via Euclidean methods; Caron-Huot showed that, at next-to-leading order in the coupling, the desired Wilson loop is equivalent to a modified Wilson loop [11] in a 3-dimensional Euclidean theory called Electrostatic Quantum Chromodynamics, or EQCD. EQCD was first postulated by Braaten and Nieto [12]. In high-temperature QCD, as described by EQCD, only the gluon Matsubara-0-mode remains dynamical, all other degrees of freedom can be integrated out and contribute to the effective field theory parameters. While the spatial components of the gluon field persist in EQCD with a modified coupling g 2 3d , the temporal field component A 0 is no longer protected by gauge invariance from acquiring a screening mass m 2 3d . Thus A 0 turns into a scalar in the adjoint representation of SU(3) with a mass m D and a quartic self-coupling λ. A perturbative matching between full QCD and EQCD has been carried out for all parameters up to at least two-loop order, for instance in [13,14].
Caron-Huot showed that, at next-to-leading order in the coupling and up to IR safe corrections, EQCD also describes correlation functions of spacelike or lightlike operators at sufficient distance, and therefore the correlation functions relevant in the Wilson loop which determines C(b ⊥ ) (at least for b ⊥ 1/2πT ). Because EQCD captures the nonperturbative IR dynamics which spoil the convergence of perturbative computations for thermodynamical quantities such as the pressure [15][16][17][18][19][20] and correlation lengths [21,22], we expect this EQCD approach to work for C(b ⊥ ) wherever dimensional reduction [23,24] works, which may be as low a twice the QCD crossover temperature [25].
We should exploit this opportunity to determine a key dynamical property of QCD, highly relevant for jet quenching [2,26] and transport coefficients [9], by performing a nonperturbative determination of C(b ⊥ ) within EQCD, on the lattice. This would be tremendously assisted by a complete determination of all O(a) corrections between lattice and continuum EQCD, both for the Lagrangian and for the Wilson loop operator in question. The renormalization up to linear order in the lattice spacing a is known analytically for almost all parameters of EQCD [27]. More recently, the renormalization of the modified, so-called 'Null Wilson lines' was determined [28]. The only missing O(a)-contribution to renormalization stemmed from the mass m 2 3d ; we recently determined it numerically in [29]. This development allows for precision studies of EQCD free from any O(a) discretization errors.
In this paper we undertake a comprehensive study of C(b ⊥ ) in EQCD on the lattice. We hope the resulting nonperturbative information about transverse momentum diffusion can be helpful in the future, to study jet energy loss and thermalization in the Quark-Gluon Plasma at temperatures where dimensional reduction is useful. The next section will define the problem more completely. Section 3 will describe in detail our lattice procedure. Section 4 will present our results, and finally we will end with a short discussion. Tabulated results and covariance matrices appear in an appendix.

Formulation of the problem
Electrostatic QCD is a dimensionally reduced effective theory for high-temperature QCD. At high temperatures T , the Euclidean path integral extends in all spatial directions but is periodic in Euclidean time with extent β = 1/T , and periodic/antiperiodic boundary conditions for integer/odd-half-integer spin fields. The long-distance physics is dominated by zero-Matsubara modes of the gauge fields, or equivalently, the most important IR physics comes from A µ fields which are constant across the temporal direction. Fluctuations which vary in the temporal direction, including all fermionic fluctuations, can be integrated out and contribute via the EFT parameters. The temporal gauge field component A 0 can be interpreted as a scalar Φ, living in the adjoint representation of SU(3). Gauge invariance no longer protects it from receiving a mass; both a mass and a quartic interaction are generated when we integrate out fluctuations, leading to a continuum action which reads The gauge coupling g 2 3d is dimensionful, and approximately equals g 2 4 T with g 2 4 the squared gauge coupling of the original 4D theory. Therefore one can consider g 2 3d to set a scale (an energy or inverse length scale), and we can use it to express the other couplings in terms of dimensionless ratios: x ≡ λ/g 2 3d , which expresses the strength of the scalar self-coupling, and y ≡ m 2 D (μ = g 2 3d )/g 4 3d , which expresses how heavy the scalar mass is. Formally x ∝ α s and y ∝ α −1 s , so when the coupling is perturbative we are in the small x and large y regime. The relation between the parameters x, y and the gauge coupling and number of light fermions of full QCD have been worked out to the two-loop level [14] and are illustrated in Figure 1. In this paper we will consider four specific cases, listed in Table  1.  Our observable of interest is the modified Wilson loop proposed in [11],   [30], and the physically relevant points lie in the supercooled "high-temperature" symmetric phase below the transition line. The blue points mark our scenarios of interest from Table 1. contain both the gauge field and the exponential of the Φ field: Note that Φ appears without a factor i, that is, the long edges of the Wilson loop are not unitary. This term can be understood as −iA 0 in the original Minkowski picture, with A 0 rewritten as Φ and rotated by a factor i when we pass to Euclidean signature. Our goal is to establish the value of C(b ⊥ ) using Eq. (2.2) for several transverse distances, in each of the four cases listed in Table 1.

Lattice implementation
The EQCD continuum action was discretized on a three-dimensional spatial lattice, including the now-complete O(a) corrections to the lattice parameters. The lattice fields are updated using a mixture of heatbath and overrelaxation updates. One full update consists of a sweep over all sites of one heatbath update to each scalar and link, followed by four sweeps in which each scalar and link is updated with over-relaxation. The presence of an adjoint scalar requires an additional accept-reject step be added to the gauge boson updates. Our code involves custom modifications of the openQCD-1.6 codebase [31]. For more details of our lattice implementation, we refer to [29].
The modified Wilson loop (2.2) on the lattice is a local observable which suffers from significantly higher noise levels than volume-averaged observables such as the scalar con- . This effect is especially strong for loops that enclose large areas, which automatically applies in our case since we would like to extract an observable that requires an extrapolation of L → ∞. An algorithm that was designed to overcome that problem is the multilevel algorithm proposed by Lüscher and Weisz in 2001 [32,33]. It relies on freezing one or multiple surfaces perpendicular to the Wilson loop and updating the subvolumes separately. This allows to average over the subvolumes independently and measuring the final observable as a correlation of multiple, pre-averaged quantities, which reduces the noise drastically. Originally, this algorithm was designed for pure SU(3) Yang-Mills theory, but it found an application to EQCD, too [34]. We split our lattices along the largest extend N z in 4 sublattices, on which we perform 80 update sweeps separately before an update sweep through the complete volume is conducted.
Following Panero, Schäfer and Rummukainen [34], the lattice implementation of the modified Wilson loop reads where U i is the standard gauge link in i-direction, the transverse separation is assumed to be n b lattice spacings a in the x-direction b = n b ae x and Z is the renormalization factor of the Null Wilson line. This quantity was analytically computed to O(a) in [28]. We repeat their main result for SU(3) explicitly for convenience: where g 2 3d a is the dimensionless lattice spacing, Z Φ is the overall normalization factor of the scalar Φ, which we set Z Φ = 1 without loss of generality, and Σ = 3.17591153562522 and ξ = 0.152859324966 are standard integrals in the lattice-continuum-matching of EQCD.
When computing a Wilson loop it is often possible to replace each link with its heatbath average in computing the Wilson loop's value. This is possible so long as no term in the lattice action contains more than one link which appears in the Wilson loop. Therefore we cannot apply this method here, because the scalar field appears in an action term containing the link variable and it also appears in the modified Wilson loop. Similarly, we have not found an efficient way to average over scalar fields, because each scalar field depends on the neighboring scalar fields and because of the quartic term in the action. Therefore the multilevel procedure is the only noise reduction technique we have applied.

Extracting
The previous section explains how we obtain data forW (L, b ⊥ ) at multiple lengths L, transverse distances b, and lattice spacings g 2 3d a, each at four "temperatures" (x, y choices). Here we explain how we use this data to extract C(b ⊥ ). There are three different limits that have to be taken into account properly: Since EQCD possesses a mass gap, the first limit can be easily reached by choosing sufficiently large volumes [35]. The continuum limit is performed by a standard polynomial extrapolation of g 2 3d a → 0, where the linear term was eliminated by the full O(a) improvement. So the remaining limit to be treated is the infinite length limit, which we will discuss in the next few paragraphs.
Unfortunately, the information about C(b ⊥ ) inW (L, b ⊥ ) is contaminated by (in principle infinitely many) higher states' energies [36] where c i are prefactors resulting from the geometry of the Wilson loop and lattice artifacts and are not important for our purpose. Fitting a large number of exponential functions is in general a very hard problem, which gets even worse if the decay constants are numerically close to each other. What saves the day is that the energies are increasingly ordered in n, ie. their exponential decay happens faster at large L as n increases and we are only interested in the lowest energy C(b ⊥ ). Conversely, the relative error of a Wilson loopW (L, b ⊥ ) scales inversely with the enclosed area, so small loops feature good statistics. It is crucial to find a regime in which the balance between sufficiently small Monte Carlo errors and sufficiently large L for small contamination is maintained. In our case, g 2 3d L ≥ 1.0 fulfilled that requirement such that it is sufficient to consider only one higher state. Thus, the fit function for the L → ∞ extrapolation reads lattice data 1−state fit 2−state fit Figure 2. Genuine L → ∞ fit for g 2 3d b ⊥ = 2.5 at g 2 3d a = 1/6 and x = 0.0677528. Note that this figure only illustrates the initial guess finding. The actual values extracted from the data by the variable projection cannot be displayed, since this procedure does not give values for the c i 's. However, we found that the final results from the variable projection were very close to the initial ones determined as in the plot.
Exponential fits of this form are in general very unstable. There are a few techniques one can apply to improve that, however. Firstly, one can estimate starting values close to the final fit values, for which a procedure was outlined in [37]. The second method we apply is the so-called variable projection [38]. Roughly speaking, one gives up on determining the c i 's and finds the minimum of χ 2 in the reduced parameter space, only. In our case, this is an appropriate procedure since we are not interested in the values of the c i 's, anyway.
As a last obstacle, we determine our data for allW (L, b ⊥ ) at a given lattice spacing from the same ensemble, which means that our data is highly correlated, not only along the Monte Carlo time axis, but also for the different L and b ⊥ . The correlation along the Monte Carlo time axis can be eliminated by binning the data; the bin size has to be varied until a plateau for the errors is reached. The correlation of different lengths is taken care of by performing correlated, variable-projected fits [38,39]. Last but not least, the correlation for the final, continuum-extrapolated different b ⊥ is less severe than the one for the different lengths since multiple lattice spacings (ensembles) contribute to the continuum-extrapolated points of C(b ⊥ ). Nevertheless, we report the covariance matrices  This data can now be extracted at various lattice spacings and extrapolated to the continuum, cf. Fig. 3. We choose the lattice spacings such that a quadratic interpolation is sufficient. Since the linear term is eliminated by the full O(a) improvement, all continuum extrapolations are fits with only two free parameters, improving the error of the extrapolated value drastically. The shape of the extrapolation fit in Fig. 3 is clearly quadratic, confirming that our improvement procedure succeeded in eliminating all linear-in-a renormalizations and rescalings.

Analytical expectations: small b ⊥
Before presenting numerical results, we should start by asking, what answers do we expect? Of course we don't know what the behavior of C(b ⊥ ) should be, otherwise there would be no need to measure it nonperturbatively on the lattice. But in limiting regimes, namely g 2 3d b ⊥ 1 and g 2 3d b ⊥ 1, we might expect simpler behavior. Let us start with g 2 3d b ⊥ 1. First note that C(b ⊥ ) has the same units as energy. To see this note that the log trace of a Wilson loop is dimensionless, so Eq. (2.2) shows that C(b ⊥ ) has units of inverse length or energy. Alternately, C(q ⊥ ) describes a probability per unit length and momentum-squared, and so has units of inverse energy. Fourier transforming d 2 q ⊥ introduces two factors of energy, making C(b ⊥ ) linear in energy. Next, note that both g 2 3d and m D have units of energy. Perturbation theory is an expansion in g 2 3d , but since this quantity is dimensionful it should be balanced by powers of something else with dimensions; hence we have an expansion in g 2 3d b ⊥ and/or in g 2 3d /m D . Therefore we formally expect that the small-b ⊥ expansion, as an expansion in loop order, is of form If we were computing the standard Wilson loop, the first two leading-order terms might arise; but for the structure we do consider, there are cancellations at leading order between A z and Φ contributions, precisely because Φ couples in an antihermitian way. Since the Φ field is heavy, the cancellation is incomplete and the g 2 3d b 2 ⊥ m 2 D term is present, but the g 2 3d and g 2 3d b ⊥ m D terms are absent. In particular, the leading short-distance contribution is [40] which as expected scales as g 2 3d m 2 D b 2 ⊥ . Note that a quadratic term in C(b ⊥ ) can be understood, using Eq. (1.1), as C(b ⊥ ) q 4 b 2 ⊥ (see [9] Appendix C). Therefore the logarithmic term here represents a log UV divergence in the Coulombic value ofq, which is well known. Also note that the dominant momentum region giving rise to Eq. (4.2) is q ⊥ ∼ m D , since this is the momentum region where the cancellation between A z and Φ first breaks down. The momentum region q ⊥ ∼ g 2 3d gives rise to an O(g 6 3d b 2 ⊥ ) contribution, which therefore indicates the order where nonperturbative physics will enter.
At NLO the full b ⊥ m D dependence has been worked out in Ref. [9], based on q ⊥ -space results from [11]. They find The linear term is the only term linear in b ⊥ which will arise, and is therefore a clean prediction of perturbation theory. The second term is formally suppressed relative to the LO expression by a factor ∼ g 2 3d /m D , indicating that in this region, perturbation theory is an expansion in (g 2 3d /m D ) ∼ y −1/2 . Similarly, the unknown NNLO contribution is of order g 6 3d b 2 ⊥ . Unfortunately a 2-loop calculation would not be sufficient to determine this term, since as we already discussed, this is the order where perturbation theory for this quantity breaks down due to IR divergences; all higher loop orders also contribute at this order, so the g 6 3d b 2 ⊥ coefficient is nonperturbative. This unknown nonperturbative entry is suppressed, with respect to the leading-order result, by a factor of y −1 . Therefore, the highest-temperature case we consider, with y = 1.65, should show reasonable convergence and the two known perturbative terms should be relatively close to determining the true linear plus quadratic behavior at small b ⊥ . But for the other values we consider, perturbative results for the b 2 ⊥ coefficient will not be useful and this coefficient (and thereforeq) has to be fitted. However the b ⊥ and b 2 ⊥ ln(b ⊥ m D ) terms are clean predictions of perturbation theory.

Analytical expectations: large b ⊥
The large b ⊥ region has been discussed by [34], who argue that g 2 3d b ⊥ 1 corresponds to the region where Wilson loops display area-law behavior. The Φ field correlator essentially vanishes between opposite edges of the Wilson line and does not contribute to the b ⊥ dependence in this regime; therefore we expect where σ EQCD ∝ g 4 3d is the EQCD string tension. The EQCD string tension was predicted in [14], continuing the continuum-extrapolated result of three-dimensional pure gauge theory [41], ie. magnetostatic QCD (MQCD) to EQCD. However, this calculation relies on a perturbative matching of the MQCD coupling to its EQCD equivalent, which is again a formal expansion in 1/ √ y which may not show good convergence. Instead, we can compute the string tension of EQCD from our simulations; by setting Z = 0 in (3.1) we recover the standard Wilson loop. The trace of a standard Wilson loop should depend on its transverse and longitudinal extents, if both are large, as where σ EQCD is the area-law or linear-confining coefficient, B is a coefficient associated with perimeter-law contributions, and A is a constant associated with the corners. Both A and B are expected to suffer from UV divergences and are therefore lattice-spacing dependent, but σ EQCD should have a valid continuum value, which can be obtained as Therefore it is straightforward to predict the large b ⊥ behavior of C(b ⊥ ), which grows linearly with a coefficient which is nonperturbative but can be independently evaluated.

Numerical results
Section 3 has already described how we extract C(L, b ⊥ ) values at each lattice spacing, and how we extrapolate these to the large L and small g 2 3d a limits. Doing so, we find the results tabulated in Table 2 and displayed in Fig. 4. These represent our principal findings.
Note that every result in Table 2 is based on data from at least three lattice spacings, extrapolated to the continuum -except for the first (smallest-separation) data point in each column, see Table 3 in the appendix. This point is based on a single lattice spacing. Analyzing the lattice-spacing dependence of the points where we can perform a continuation, we find that the a 2 coefficient is fairly constant for small b ⊥ ; so we assume that the same a 2 extrapolation coefficient applies for this smallest-b ⊥ point as for the next-smallest b ⊥ data, and assign 100% systematic errors to this estimate. We then combine this (systematic) error with the statistical error in quadrature.
x cont = 0.08896 x cont = 0.0677528 x cont = 0.0463597 x cont = 0.0178626 y cont = 0.452423 y cont = 0.586204 y cont = 0.823449 y cont = 1.64668  Table 2. Results for C(b ⊥ ) for four temperatures and a range of transverse separations. All data points are continuum extrapolated, with errors representing all statistical and systematic errors associated with the data extraction and extrapolations, except for the first (smallest g 2 3d b ⊥ ) entry in each column; see text. We also quote the extracted value ofq and the string tension as determined from Eq. (4.6), see text.
Let us examine the small b ⊥ region in more detail. As we saw in Subsection 4.1, we expect that C(b ⊥ ) scales, for small b ⊥ , as (4.7) Figure 5 plots the small b ⊥ data with two fits. The long-dashed curves are based on assuming thatq is determined by the NLO expression presented in Eq. (4.3) -that is, neglecting the (nonperturbative) g 6 3d corrections. The dashed curves represent a fit of the data, using all data points with g 2 3d b ⊥ ≤ 0.75, to the functional form shown in Eq. (4.7), treatingq as a free fitting coefficient and neglecting O(b 3 ⊥ ) effects (a one-parameter fit). The resulting value forq (which, note, is corrected by the ln(b ⊥ m D ) term), appears as an added line of Table 2. This result is about a factor of 2 smaller than the NLO result, indicating that the next correction is never small.
Finally, we consider the large b ⊥ asymptotics. As discussed, one expects C(b ⊥ ) at large b ⊥ to rise linearly, with a coefficient which equals the string tension. We determine the string tension numerically using the procedure outlined in Subsection 4.2, specifically Eq. (4.6), for each (x, y) pair (temperature), reporting our results in the last line of Table  2. These determined slopes are then compared to the C(b ⊥ ) results in Figure 6. The figure shows that the large-b ⊥ asymptotics are indeed well described by the string tension in EQCD. The range over which this linear behavior holds is larger for the high-temperature (small-x, large-y) case, where the scalar decouples over a shorter distance. For our highest-  temperature case we find nearly linear behavior already close to g 2 3d b ⊥ 1. For the lowest temperature, on the other hand, the scalar is very light and its effects persist to larger distances; here it takes until approximately g 2 3d b ⊥ 4 before the asymptotic behavior sets in.

Conclusion and outlook
The collision kernel C(b ⊥ ) contains essential information about how a thermal medium transverse-broadens, and therefore damps, high energy particles. Its perturbatively illbehaved infrared contributions can be computed within the effective theory EQCD by formulating the collision kernel in terms of EQCD variables [11] and computing them on the lattice [34]. We have presented the first such lattice analysis of C(b ⊥ ) which is complete in the sense that it uses all lattice-continuum improvements and makes a complete extrapolation to the continuum limit. Our approach profited from the use of the multilevel algorithm for noise reduction and the variable projection method for fitting to the long Wilson-loop limit. The use of fully improved lattice-continuum matching and operator definitions accelerated the continuum limit, leading to high precision continuum-extrapolated results.
Our results indicate a rather large downward correction in the value ofq relative to the NLO result; indeed, we find aq value which is closer to the leading-order result. This  has implications for jet quenching and for transport coefficients, which are also sensitive toq. However we want to emphasize here that the reader should not dwell on the value ofq itself, which is a crude fit to a few small-b ⊥ points. Rather, interested practitioners should directly use C(b ⊥ ), which is the potential, in the transverse plane, which is relevant for driving decoherence within Zakharov's spacetime picture of jet quenching [42][43][44][45].
Another interesting feature of our result is that the value of the kernel C(b ⊥ ) is actually negative for sufficiently small g 2 3d b ⊥ values. This behavior is both expected and confusing. In Ref. [9] it is shown that C(b ⊥ ) should have a small-b ⊥ expansion displaying a leading negative, linear-in-b ⊥ behavior. A negative value of C(b ⊥ ) does not have a sensible probabilistic interpretation; but this negative linear term can only dominate the result for b ⊥ ∼ g 2 3d /m 2 D , which is formally O(1/T ). This is exactly the short-distance regime where EQCD breaks down as an effective description of thermal QCD. Unfortunately, short distance, which corresponds to high transverse momentum, is relevant for the highest-energy splitting processes, which may be physically important in the medium-modification of high energy jets. Therefore it appears to us that more effort needs to be focused on obtaining a better understanding of the matching the behavior of C(b ⊥ ), especially at short distances, between EQCD and full QCD, so that our results can be correctly incorporated into jet-quenching calculations. Such an analysis would also shed more light into the role of collinear effects, which are expected to enter at the NNLO level but should be enhanced by double logarithms [46] and may involve more complex structures than the Wilson loop considered here [47,48]. We leave this, and an application of our results to the computation of jet modification, for future work.

A Simulation parameters and correlation matrices
For completeness, we provide a "data dump" of the details of our simulated boxes, statistics, and correlation matrices. Because the correlation matrices are symmetric and nearly banddiagonal, we will only list the entries 1 to 3 places to the right of the diagonal.
We also explain what at first sight is an odd choice of lattice volumes. The volume should be chosen such that N x/y/z ≥ β in order to avoid finite volume effects [35]. Furthermore, it has to fulfill N x/y > 2b max ⊥ and N z > 2L max to suppress an (unphysical) correlation over the boundaries of the simulated box. Since not all g 2 3d a lattices produce sensible information at all g 2 3d b ⊥ and only 3 lattices are required for a continuum limit, the choice of volumes seems a little odd at first sight.      Table 7. Correlation matrix for x = 0.0178626 case.