Probing the photon emissivity of the quark-gluon plasma without an inverse problem in lattice QCD

The thermal photon emissivity of the quark-gluon plasma is determined by the in-medium spectral function of the electromagnetic current at lightlike kinematics, σ ( ω ). In this work, we present the first lattice QCD results on moments of σ ( ω ) /ω , defined by the weight function 1 / ( ω 2 + (2 πT n ) 2 ), n ∈ Z and computed without encountering an inverse problem. We employ two dynamical flavours of O( a )-improved Wilson fermions at a temperature T ≈ 250 MeV and perform the continuum limit. We compare our results for the first two moments to those obtained dispersively by integrating over the spectral function computed at weak coupling by Arnold, Moore and Yaffe.


I. INTRODUCTION
The strongly-interacting matter created during a heavy-ion collision radiates electromagnetic energy due to the interactions of its charged constituents.Photons, the quanta of this radiation, have a large mean free path compared to the strongly-interacting degrees of freedom, thanks to the hierarchy of the electromagnetic and strong coupling constants at experimentally accessible temperatures, α em ≃ 1/137 ≪ α s .These probes escape from the strongly interacting medium and can inform us about its entire space-time evolution.
Let ω be the energy of a photon emitted from a fluid cell at rest and in thermal equilibrium.The differential rate of photon emission per unit volume of the cell is described to all orders in the strong coupling, and to leading order in the electromagnetic coupling by the formula [1] where σ(ω) is the in-medium spectral function associated with the two-point correlation function of the electromagnetic current at zero-virtuality.Therefore, this temperature-dependent spectral function directly determines the emissivity of real photons from thermal matter described by quantum chromodynamics (QCD).From the experimentally measured spectrum of photons, one would hope to infer interesting physical properties of the thermal medium, starting with its temperature.However, since the QCD matter created in a heavy-ion collision does not remain in equilibrium, the prediction of the total yield is complicated by the necessity to integrate over the spacetime history of the expanding fluid as well as to include non-thermal components [2].Nevertheless, recent experimental results from the RHIC and the LHC facilities have been compared to phenomenological models which incorporate a thermal component produced by the expanding plasma [3].Data for the yield of direct photons from the PHENIX and ALICE experiments at RHIC and the LHC, respectively, show an excess of production in the region of phase space where thermal sources dominate, i.e. for transverse momenta, p T , below a few GeV [3][4][5][6].In contrast, data from the STAR experiment at RHIC appear to be in agreement with the available theoretical models [7].Furthermore, the anisotropy of the photon momentum distribution, the photon v 2 , is larger for all experiments than models suggest [8][9][10].
The failure of most models in describing collectively the experimental findings mentioned above is usually referred to as the direct photon puzzle [2,3,11].The model calculations rely on theoretical predictions of the thermal photon emissivity, which is directly proportional to σ(ω).So far, the phenomenological models have used such predictions obtained at leading order in QCD perturbation theory [12,13], complemented with relativistic kinetic theory calculations in a hot meson gas for the hadronic phase (see e.g.[14,15]).How exactly to connect these predictions in the vicinity of the crossover temperature of about 150 MeV is unclear, hence interpolations have been used.
Lattice QCD is a suitable framework for the non-perturbative determination of the correlation function which mathematically determines σ(ω) uniquely.However, given that lattice QCD is formulated in Euclidean space, obtaining the spectral density by analytic continuation from the current correlator is a numerically ill-posed problem (see for instance the reviews [16,17]).The analytic continuation needed to determine real-time quantities from the correlators computed in Euclidean spacetime is further hampered by the fact that a relatively small number of noisy data points are available, usually O(10 − 50) at finite temperature.Yet what one is interested in is determining a spectral function at O(1000) values of ω covering the interesting kinematical region with a sufficiently high resolution and with a few-percent precision.
Aiming at a goal which is within reach of the current numerical capabilities, one can try to extract not the spectral function itself, but a smeared version thereof.In the method originally proposed by G. E. Backus and J. F. Gilbert [18,19], the calculation of a smeared, filtered spectral function is viable after determining certain coefficients that multiply the Euclidean correlator data points.The value at a given ω of this smeared spectral function is a weighted average of the actual spectral function over the vicinity of ω.The weighted average is defined as an integral over the spectral function multiplied by a compact kernel centered at a certain value of ω.This method has been applied to several problems [20][21][22][23][24][25][26][27] in lattice QCD, and similar approaches have been proposed which offer a path to reducing the width of the filtering kernel, thus getting closer to the unsmeared spectral function [28][29][30].These approaches have in common that they provide information about an integral over the spectral function based on the Euclidean correlator and additional technical ingredients.
In this work, we obtain integrals over σ(ω) by calculating a correlator directly accessible in lattice QCD that has be shown [31] to have a simple integral representation in terms of σ(ω).The key ingredients are to write a dispersion relation for the Euclidean correlator at fixed photon virtuality, rather than at fixed spatial momentum, and to employ an imaginary spatial momentum in order to realize lightlike kinematics in Euclidean space.We thus obtain information about the spectral function at lightlike kinematics without facing an inverse problem.This cannot be achieved when employing dispersion relations at fixed spatial photon momentum, which has been the only type of dispersion relation applied to date in lattice-QCD based studies of non-equilibrium properties of the quark-gluon plasma, starting with Ref. [32].
Specifically, we compute the first two moments of σ(ω)/ω, defined by the weight function 1/(ω 2 + (2πT n) 2 ) for n = 1 and n = 2.The construction makes use of the spatially transverse Euclidean correlator evaluated at Matsubara frequency ω n and at imaginary spatial momentum, k = iω n .Such an 'imaginary momentum' corresponds to a weight function exp(ω n x 3 ) in the spatial coordinate x 3 , which strongly enhances the contribution from large positive x 3 .Thus the control over the statistical fluctuations at long distances is the main difficulty in computing these observables in lattice QCD.Numerically, the task thus bears a strong resemblance with the lattice calculation of the hadronic vacuum polarization contribution to the muon g − 2 [33][34][35].Using the spectral function σ(ω) computed in QCD at weak coupling [12,13], as well as that computed in strongly coupled super-Yang-Mills by AdS/CFT methods [36], we compute the same moments and compare them to our results.It is the first time that the weak-coupling QCD predictions for the photon emissivity of the quark-gluon plasma can be tested non-perturbatively without any uncertainties associated with an ill-posed inverse problem.
The structure of the paper is the following: In Sec.II A, we introduce the basic observables in the continuum.In Sec.II B, we discuss various subtractions that are necessary to reduce large cutoff effects afflicting our primary observable, while Sec.II C discusses the virtuality dependence of this observable.The lattice observables are defined in Sec.II D. Thereafter, in Sec.II E, we provide some details about our numerical setup.This is followed by the presentation of the main results in Sec.III and Sec.IV, in which we devote different subsections to the different derived observables that we investigated and describe our analysis for the sectors n = 1 and n = 2, respectively.In Sec.V, we compare our results to predictions made using either free quarks, the weak-coupling expansion or the strongly-coupled N = 4 super Yang-Mills theory (SYM).Finally, we conclude in Sec.VI.

II. PRELIMINARIES A. Definitions in the continuum
The spectral function of the electromagnetic current is defined as Here, J em µ (x) = f Q f ψf (x)γ µ ψ f (x) is the electromagnetic current, Q f ∈ {2/3, −1/3} denotes the charge of a quark with flavor f and γ µ are the Euclidean Dirac matrices satisfying {γ µ , γ ν } = 2δ µν , i.e. γ µ † = γ µ .The time evolution is given in real time by J em µ (x) = e iHt J em µ (0)e −iHt .The expectation value is taken with respect to the thermal density matrix e −βH /Z, where Z is the grand canonical partition function, H is the QCD Hamiltonian and β = 1/T is the inverse temperature.The thermal photon production rate per unit volume of the QGP, dΓ γ (ω)/dω, is related to the transverse channel spectral function evaluated at lightlike kinematics, σ(ω) ≡ ρ T (ω, k = ω), where is the transverse channel spectral function.
We remind the reader that Eq. ( 1) is valid at leading order in the electromagnetic coupling constant, but to all orders in the strong coupling constant [1].
Let us consider now the Euclidean screening correlators1 with ω n = 2nπT being the nth Matsubara frequency and the time-evolution being given in Euclidean space-time (O(x 0 ) = e Hx 0 O e −Hx 0 ).In the presence of interactions, the lowlying spectrum of the spectral function corresponding to Eq. ( 4) is gapped and discrete, resulting in an exponential falloff for G E, µν (ω n , p 1 , p 2 ; x 3 ).The corresponding energies are called screening masses and have been investigated in Ref. [37] using weak-coupling theory as well as lattice simulations.We discuss these in Sec.III E and in App.F.
Restricting the further discussion to the transverse channel, we define and two special cases of Eq. ( 5), the non-static (ns) and static (st) screening correlators in the transverse channel as and respectively.Thus, for the non-static (static) screening correlator, the momentum is inserted into the temporal (spatial) direction.The Fourier transform of the non-static screening correlator is defined as Evaluating Eq. ( 8) at imaginary spatial momentum, we define [31] the spatially transverse Euclidean correlator evaluated at Matsubara frequency ω n and at imaginary spatial momentum k = iω n .Correspondingly, we can also write An important property of H E (ω n ) is that it vanishes identically in the vacuum, a consequence of Lorentz symmetry [31].Secondly, subtracting explicitly H E (0), which vanishes due to current conservation, suffices to make the expression on the right-hand side of Eq. ( 11) finite by power counting.Indeed, this subtraction amounts to the modification e ωn(ix 0 −x 3 ) → (e ωn(ix 0 −x 3 ) −1) inside the integral, and Taylor-expanding the exponential function yields the leading term 1 2 ω2 n (x 2 3 − x 2 0 ) if we drop terms that do not contribute due to the x 3 → −x 3 or x 0 → −x 0 symmetry of the ⟨J 1 (x)J 1 (0)⟩ correlator.The terms 1  2 ω 2 n (x 2 3 −x 2 0 ), which separately would generate a logarithmic UV-divergence, obviously cancel each other exactly in the vacuum, and therefore only make a UV-finite contribution at non-zero temperature.This last statement is easily proven by power counting, using the operator-product expansion.Even though the simple subtraction (e ωn(ix 0 −x 3 ) −1) is thus in principle sufficient, in subsection II B we will derive other possible subtractions, that turn out to be superior in practical lattice QCD calculations, and test them in the context of lattice perturbation theory (appendix E), In Ref. [31] it has been shown that the spectral function at vanishing virtuality, σ(ω), is related to H E (ω n ) via the following once-subtracted dispersion relation The left-hand side contains the difference of two observables that can be evaluated in Euclidean space-time and this difference directly probes an integral of the spectral function at vanishing virtuality.
Given that H E (0) = 0 in the absence of massless static screening modes, a simplified dispersion relation applies for A question arises as to how to interpret the right-hand side of this equation for ω n = 0.The rule that applies to all cases of interest in the quark-gluon plasma context is that one should simply interpret it as being zero for ω n = 0.In an interacting theory 2 , where we expect σ(ω)/ω to be finite at ω → 0, the right-hand side also vanishes if one takes its limit for ω n → 0; however, in the theory of free thermal quarks where σ(ω)/ω ∝ δ(ω), the right-hand side would not vanish in that limit.This issue amounts to the same question as to whether the function H E (ω n ), analytically continued to all frequencies, is continuous at argument zero (see [16], sec.2.3).The integration kernels multiplying the function σ(ω)/ω in Eqs. ( 12) and ( 13) are shown in Fig. 1.One can notice the different characteristics of these curves, the integrand of Eq. ( 13) being sensitive to soft photon emission, while the kernel for H E (ω 2 ) − H E (ω 1 ) tends to zero for ω/T → 0. This difference is significant, since the slope of σ(ω) at the origin gives access to the charge diffusion coefficient D [23,38,39],  12) and ( 13) multiplying the function σ(ω)/ω.Right: The spectral function normalized by 2χ s ω/T at lightlike kinematics in QCD at complete leading order according to Ref. [12,13], and in the strongly-coupled N = 4 SYM theory [40]. where is the static charge susceptibility.We recall that for free massless quarks,

B. Lattice subtractions
The lattice regulator breaks Lorentz symmetry: short-distance contributions to H E (ω n ) emerge which would hamper the continuum extrapolation.It is therefore necessary to develop suitable lattice representations of H E (ω n ).In fact, for discretizations of the correlator which use two local currents or two exactly-conserved lattice currents, such subtractions are required to cancel divergences arising from contact terms and have a well-defined continuum limit, see Sec.II D for further discussion.This can be achieved in several ways, for instance by subtracting the vacuum lattice correlator obtained at the same bare parameters, as was proposed in Ref. [31] or by subtracting a thermal correlator having the same short-distance properties [41].
The basic observation is that the Fourier-transform of a static screening correlator at lightlike momentum vanishes in the continuum.Indeed, the restriction of the polarisation tensor to spatial components in the static sector has the form familiar from the vacuum polarisation.From here and from the absence of a pole in Π at p 2 = 0 follows immediately the property Thus, generalizing the estimator proposed in [41], we consider a subtraction involving the static screening correlator at momentum p, We note again that in the case of the static (st) transverse channel screening correlator (defined in Eq. ( 7)), the momentum p is inserted into a spatial direction (here x 2 ) orthogonal to x 3 and to the directions corresponding to the Lorentz indices of the currents (here x 1 ).The case of p = ω n [41] has the special property that H E, ωn (ω n ) vanishes identically at T = 0, correctly reproducing the continuum property that H E (ω n ) vanishes in the vacuum.This property is expected to reduce discretisation errors of H E, ωn (ω n ) at non-zero temperature, an expectation that is confirmed in lattice perturbation theory (see App. E).
We remark that even more general subtractions are possible, i.e. one can subtract a general linear combination of static screening correlators at different values of p i .Eq. ( 18) is a special case of Eq. ( 19) with α 1 = 1, p 1 = p and α i>1 = 0.When evaluated at different values of p, the results for H E, p (ω n ) in Eq. ( 18) do not agree with each other in general on a finite lattice, but they have to match after taking the continuum limit.The same holds for the quantity in Eq. ( 19), i.e. choosing different α i coefficients and subtracting the static screening correlators evaluated at different momenta, the results for H E, p, α (ω n ) differ at a finite lattice spacing, but have to agree in the continuum.One can exploit these observations and propose subtractions that may have more tractable integrands and/or reduced cutoff effects toward the continuum.Similarly, a direct probe of the difference of H E (ω n ) and H E (ω r ) on the lattice is provided by the discretized version of Again, the contribution coming from the static screening correlators in Eq. ( 20) vanishes in the continuum.We will explicitly investigate the more general subtracions for the extraction of H E (ω 2 ), i.e. in the second Matsubara sector, as well as for the difference H E (ω 2 ) − H E (ω 1 ) in Sec.IV.Additionally, alternative subtractions concerning H E (ω 1 ) and H E (ω 2 ) will be exploited in Apps.B and C, respectively.In the rest of the paper we use the notation H E omitting the extra subscripts to indicate the standard subtraction with α 1 = 1, α i>1 = 0 and p = p 1 = ω n .

C. Probing the virtuality dependence of H E
Until this point, the above discussion focused on observables at vanishing virtuality, since this is the relevant kinematics for the photon emissivity.We note that H E can also be evaluated for a given non-vanishing virtuality, which can be helpful to understand the behavior of the low-mass dilepton rate [42].In order to investigate the effect of introducing a small non-vanishing Q 2 , we can calculate the derivative of H E with respect to Q 2 , evaluated at Q 2 = 0.For that, we start from the definition given in Ref. [42] By expanding around Q 2 = 0 for ω n > 0 we obtain Evaluating the derivative in the free theory, we can use the asymptotic formula G T E,ns (ω n , x 3 ) x 3 →∞ ∼ e −ωn|x 3 | /x 2 3 (see Ref. [37], Eq. (3.11)).With this we find that in the free theory, the O(Q 2 ) term contains an infrared-divergent term, with c 1 = N c /(32π 2 ).This coefficient can be obtained by introducing a cutoff 1/Q for the x 3 -integral.
In order to ensure that we have a definition that is ultraviolet-finite, we remove the divergence present at finite Q 2 by subtracting H E (0, Q 2 ), which does not change the value at Q 2 = 0 of the function we are Taylor-expanding because H E (0, 0) = 0. Thus we evaluate In Eq. ( 24), G T st (0, x 3 ), which denotes the static, zero-momentum screening correlator, does not yield an infrared-enhanced contribution.

D. Lattice observables
In this Section, we introduce the lattice observables that we have investigated.We define the bare local and the conserved vector current as and (26) respectively, where Ψ = (u, d) ⊤ represents the isospin doublet of mass-degenerate quark fields and τ 3 is the diagonal Pauli matrix.In other words, we focus on the isovector flavour combination, the current being normalized according to f Q 2 f = 1.It is in that normalization that our results for H E (ω n ) are given.For an estimate of the physical photon emissivity assuming SU(3) flavour symmetry among the (u, d, s) quarks, one must include the factor 25) and ( 26), we define the following bare, non-static screening correlators: The Greek letters α and β stand for the discretization of the current at sink and source, respectively.We call the above correlators non-static, because the momentum is inserted into the Euclidean time-direction.We denote these with the subscript ns.By contrast, when injecting the momentum into a spatial direction, k, -perpendicular to direction of x i -, we got the bare, static screening correlators at finite momentum, (28) After choosing a particular spatial decay direction (direction of the correlator separation), i, we obtain the screening correlators in the transverse channel by choosing µ = ν orthogonal to this decay direction.We note here that we do not discriminate the notation used for the continuum or the lattice observables, i.e. we use capital G for the lattice and for the continuum screening correlators as well.
We also specify the correlator at a momentum transverse to both µ = ν and i.Therefore, we have in total six possible combinations of the decay direction and of the Lorentz indices of the currents for the non-static screening correlator and also six combinations of the decay direction, the Lorentz indices of the currents and the momentum inserted for the static screening correlator.We average over these different screening correlators measured on the same configuration.Moreover, the local-conserved and conserved-local discretizations can be transformed into each other, using Cartesian coordinate reflections.Therefore, we average these two and refer to this averaged correlator with the superscript LC in the following.We renormalize the correlators by multiplying by Z V (g 2 0 ) whenever the local vector current is included using the vector current renormalization constant from Ref. [43].Instead of the electromagnetic current, we use the isovector vector current whereby disconnected contributions are absent.
As mentioned already at the beginning of Sec.II B, one has to pay special attention when formulating a lattice estimator for H E , since a naive implementation of Eq. ( 11) on the lattice, could lead to ultraviolet divergences.A crucial point at small separation x is the removal of the x-independent part multiplying the expectation value of the product of currents in Eq. (11).The simplest way of doing this is to subtract the static screening correlator at vanishing momentum, or to subtract it evaluated at the same finite momentum ω n ,  I: Overview of the N f = 2 ensembles used in this study.Simulations were carried out at a fixed temperature of T ≈ 254 MeV and fixed aspect ratio L s /L t = 4, where L s (L t ) is the spatial (temporal) linear size of the lattice.The parameters given are the bare gauge coupling g 0 , the Wilson hopping parameter κ, the temporal size in units of the lattice spacing a, the number of configurations used N conf and the number of molecular-dynamics time units (MDUs) separating these configurations.The number of point sources per configuration is 64 in all cases.
where we used the trapezoid formula when discretizing Eq. ( 18) with p = ω n and L t (L s ) stands for the temporal (spatial) size of the lattice.The lattice spacing is deenoted by a and h(ω n , x 3 ) is the integrand, We note that G T ns as well as G T st are negative in our conventions and since in absolute value we found G T ns larger than G T st , H (sub) E (ω n )/T 2 will also be negative.From now on we leave the (sub) superscript from H (sub) E (ω n ), and similarly to the correlators, we do not discriminate between the continuum and lattice observables, i.e. we use the same symbol.
Analogously as in the case of H E (ω n ), one can obtain the lattice formula for the derivative of H E (ω n , Q 2 ) w.r.t.Q 2 by applying the trapezoid formula to Eq. ( 24), but where h(ω n , x 3 ) of the right-hand side of Eq. ( 29) is replaced by We note that this subtraction involves the completely static -zero-momentum -screening correlator, G T st (p = 0, x 3 ).

E. Simulation details
To calculate the screening correlators which enter the expression (29) for H E , we used three ensembles generated at the same temperature T ∼ 250 MeV in the high-temperature phase.We employ two-flavor O(a)-improved dynamical Wilson fermions and the plaquette gauge action; further details regarding the lattice action we used can be found in Ref. [44].The configurations for the W7 ensemble have been generated using the openQCD-1.6package and the ones of O7 and most of those of X7 using openQCD-2.0[45].512 configurations of the X7 ensemble were generated using the MP-HMC algorithm [46] in the implementation described in Ref. [47].The pion mass in the vacuum is around m π ≈ 270 MeV [44,48], and the lattice spacings are in the range of 0.033-0.05fm [44,49].The boundary conditions are periodic in space, while those in the time direction are periodic for the gauge field and antiperiodic for the quark fields, as required by the Matsubara formalism.The integrand (Eq.( 30)) for the calculation of H E (ω n ) on our finest ensemble X7, using the conserved-conserved discretization of the current-current correlator.

A. The integrand for obtaining H E
As we have shown in Sec.II D, the crucial ingredients for calculating H E through Eq. ( 29) are the non-static and static transverse screening correlators at finite spatial momentum.These have been first investigated in detail in weak-coupling theory complemented with lattice QCD simulations on a single ensemble in Ref. [37].In that work, however, the static screening correlators have been studied only at vanishing momentum.
For the determination of H E , we first analyzed the screening correlators measured on the three ensembles.We discarded a few outliers from the dataset (see App. A) and estimated the statistical errors using jackknife resampling with 200 jackknife samples.Then we formed the integrand (see Eq. ( 30)), which is shown in Fig. 2 for our finest ensemble.
As Fig. 2 shows, the quantity H E (ω 1 ) receives dominant contributions from the xT ≲ 1 − 1.5 region and with the present statistics, we have good control over the signal.The relative error on the integrand for H E (ω 1 ) is below 1% for xT ≲ 1.1 on the finest ensemble shown in Fig. 2.However, trying to evaluate H E in the n = 2 Matsubara sector reveals that one faces a severe, exponential signal-to-noise problem.The variances on the non-static correlator values are an order of magnitude larger in the n = 2 sector than for n = 1, and multiplying by cosh(2πT n) makes the situation even worse.Therefore, we focus first on evaluating H E in the n = 1 Matsubara sector and return to n = 2 only in Secs.IV A and IV B.

B. Modelling the tail of the screening correlators
While the dominant contribution to H E (ω 1 ) comes from short distances, the long-distance contribution is also non-negligible.However, it is noisier, because the screening correlators are less precise at large distances and the difference is less smooth than for short distances.Directly performing the sum of Eq. ( 29) for H E (ω 1 ) results in having relative errors of about 1.37, 1.50, 1.90% for the LL, LC and CC discretizations, respectively, on our coarsest and also 'noisiest' ensemble.We aim at a more precise determination of H E (ω 1 ) and as we will see in this section, by proper handling of the tail these errors could be reduced to 0.91, 0.77, 1.03%, respectively.
Moreover, we have an exponentially growing weight function multiplying the difference of the correlators and this results in a small enhancement of the integrand in the region x 3 ∼ L/2.This effect is the consequence of calculating the integrand in a finite volume.In order to correct for it, we applied a simple model based on fits on the screening correlators to describe the tail of the integrand.This way we have better control over the long-distance contribution and also we could correct for finite volume effects.We use the fact, that the non-static screening correlators have a representation in terms of energies and amplitudes of screening states in the following form [31]: A similar expression holds for the static correlator: The low-lying screening spectrum can be studied using weak-coupling methods as well [37].
The lowest energy of a screening state in a given Matsubara sector with frequency ω r is often called the screening mass and is denoted by In this section, we consider only the first Matsubara sector, ω 1 = 2πT , therefore in the following we do not write out explicitly the momentum dependence.
Using the above formulae, our procedure to get a better handle over the integrand is the following: 1. We split the integrand, h(x) ≡ h(ω 1 , x) (c.f.Eq. ( 30)) into two parts using a smooth step function, We call the first term in Eq. ( 34) the short-distance and the second term the longdistance part of the integrand.
2. We integrate the short-distance contribution using the trapezoidal formula.
3. We perform single-state fits on the tails of the screening correlators using the representations given in Eqs. ( 32) and ( 33) translated to a form corresponding to a periodic lattice, namely and for the non-static and for the static screening correlators, respectively.In Eqs. ( 36) and (37), L s is the spatial length of the lattice.) and G T st (ω 1 , x 3 ) respectively, using a fit range of ten lattice spacings.The bands represent the 'cosh' effective masses.Right: the tail of the integrand needed for the calculation of H E at the first Matsubara frequency.The red points are the actual datapoints on our second coarsest lattice and the band shows the result of the modelling.

Using the fit results
A ns,0 and E ns,0 as well as A st,0 and E st,0 , we replace the longdistance part of Eq. ( 34) by the corresponding infinite volume formula: x + e −(E st,0 −ω 1 )x Θ(x, x w , ∆). ( This we can integrate analytically using: where 2 F 1 denotes the hypergeometric function. While the single-state fits describe the actual data well, i.e. with good χ 2 and p-values, we note that the identification of the plateau region was not clear in some cases, although we performed a thorough scan using all possible fit ranges having different starting points and different lengths with 6a-11a.Therefore, we also made an attempt to fit the data with two-state fits with or without using priors from weak-coupling theory (c.f.Sec.F), but these fit results were not satisfactory.Typically on the coarsest or the coarser two ensembles, they either failed to describe the data, gave too large errors or returned a near-zero or negative |A| 2 coefficient -which we did not constrain -for the excited state.When the two-state fits were able to describe the data well, the ground-state static screening energy they returned was too small -as we could deduce it using the zero-momentum correlators.Therefore we decided to stick to single-state fits.Right: AIC-weighted histogram of the continuum extrapolated results for −H E (ω 1 )/T 2 obtained with x w T = 1.1.The AIC-weighted histogram of the long-distance contribution to H E (ω 1 ), shifted with the continuum result for the short-distance contribution, is also shown for comparison (sd+ld).The dark-grey band shows the systematic error, and the lighter gray band represents the total error obtained from the statistical and systematic errors added in quadrature.
Besides fitting, we also determined the effective mass using two consecutive correlator datapoints, by solving the algebraic equation for m eff .Here, G(x) denotes the actual lattice data for the non-static or the static screening correlators.We found that the effective masses are in quite good agreement with the fitted masses, but also do not show a clear and long plateau as x increases, see Fig. 3, left panel.Therefore, instead of fitting a constant, we decided to choose three representatives from a histogram built by assigning Akaike-weights [35,50] to all the fitted masses that we obtained in the most plateau-like region.In all cases, we chose a wide region, before the noise gets too large on the fitted masses.For instance for the W7 ensemble, we choose fit results on conserved-conserved correlator data in the range from x/a = 21 to 25 for both correlators (left panel of Fig. 3).We propagate the median as well as the values near the 16th and 84th percentiles of the histograms to the later steps of the analysis of H E (ω 1 ).

C. Continuum extrapolation of H
As we described in Sec.III B, we used single-state fits to describe the tail of the non-static and static screening correlators and for each ensemble and discretization we built a histogram of the fit results from the plateau region, from which we chose three representatives.When proceeding this way for each correlators, we obtain 3 × 3 = 9 possibilities for modelling the tail of the integrand, Eq. ( 38), for a given ensemble and a given discretization.We calculated H E using all these nine combinations for the tail, sorted the results and then  chose the median, the values near the 16th and near the 84th percentile after assigning uniform weights for these slightly different values of H E .This way we had three representative values of H E for each ensemble and discretization that went into the next step of the analysis, which was the continuum extrapolation.We used these in all possible combinations when performing a correlated simultaneous continuum extrapolation using a linear ansatz in a 2 .These gave (3 3 ) 3 = 19683 different continuum extrapolations.We make further variations by omitting one of the coarsest datapoints from the extrapolation, which also lead to 3 We then built an AIC-weighted histogram from using all 2 • 19683 continuum extrapolations to estimate the systematic error.A representative continuum extrapolation as well as the AIC-weighted histogram are shown in Fig. 4, left and right panel, respectively.
The transition to the modelled tail has been introduced smoothly by using a smooth step function of Eq. ( 35) and we investigated the effect of choosing different switching points, x w , in the range x w T =0.9-1.3.We found that the results were stable against these choices, see Fig. 5.
Our final result for H E in the first Matsubara sector is In order to retrieve information about the Q 2 -dependence of our observable, we evaluate the Q 2 -derivative of the difference, The continuum observable is defined in Eq. (24), and the corresponding integrand on the lattice is introduced in Eq. (31).It is interesting to have a look on this integrand -left panel of Fig. 6 -which is more pronounced at short distances than the integrand h(ω 1 , x), that we had shown for H E (ω 1 ) in Fig. 2. It starts from zero quadratically in x, and after having a peak, it crosses zero around xT ∼ 1, but the long-distance contribution is much more suppressed than it was for H E (ω 1 ).
, defined in Eq. ( 31) for the derivative of For the continuum extrapolation we applied a similar proceduce as for We remark that the result is on the order of −N c /(2π) 2 and does not exhibit any strong infrared enhancement as would be expected at very weak coupling (see Eq. ( 23)).

E. Continuum extrapolation of the screening masses
As we have already mentioned in Sec.II A, the screening masses extracted from the correlators we investigate here can also be determined using weak-coupling theory.The first relevant study in this direction was Ref. [37], which also compared the leading-order (LO) and next-to-leading order (NLO) screening masses to lattice results.However, Ref. [37] investigated only screening correlators at finite temporal or at vanishing momentum. 3Moreover, the lattice investigation has used only coarse ensembles without taking the continuum limit.
In this work, we extend this first investigation in several aspects: we investigate also screening masses obtained from screening correlators at finite spatial momentum besides the ones obtained at finite temporal or at vanishing momentum.We accumulated much more configurations and performed more measurements, enabling us to improve the signalto-noise ratio at the tails of the correlators.Finally, we also extrapolate our results to the continuum using three ensembles, from which the finest has a lattice spacing about 2/3 the lattice spacing of Ref. [37].
We start by discussing the zero and finite spatial momentum results.We performed single-state fits -discussed also in Sec.III B -using the fit ansatz given in Eq. (37).The screening correlators at vanishing momentum are more precise above xT ∼ 1.2 than at the first spatial Matsubara momentum, making it possible to determine the screening masses at 0.15-0.45%precision depending on the ensembles and discretizations.By using the dispersion relation, we could compare the extracted masses at vanishing momentum to those at momentum 2πT .With this comparison, we observe that the screening masses determinded via the dispersion relation, m 2 , give a slightly smaller value than the corresponding result at the first Matsubara momentum.
Since the single-state fitting procedure starting at later and later datapoints provides fit results converging to an asymptotic value more reliably than in the cases with finite momentum and the errors are comparable in the plateau region, we chose only one representative for each ensemble and discretization and performed the continuum extrapolation linear in a 2 using that.We also assigned a systematic error to the continuum extrapolation by removing one of the discretizations on the coarsest ensemble.The continuum extrapolation using all datapoints is shown in Fig. 7. Our continuum estimate having about 0.37% error is Using this result and the dispersion relation, an estimate for static screening mass at (spatial) momentum ω 1 = 2πT is The static screening masses being directly available at this momentum, we can extrapolate those to the continuum and see how close we get to the estimate given in Eq. (44).Therefore, we calculated a weighted average of the masses obtained by fitting the static screening correlator.We performed the averaging using the difference of smooth step functionsintroduced in Eq. ( 35) -with parameters x w T = 0.975 and x w T = 1.125, i.e. the averaging window has a width about 0.15 and was centered at 1.05.T ∆ was chosen to be 0.01.We made similar variations as before, leaving out one datapoint from the coarsest ensemble and arrived at the continuum estimate: The continuum extrapolation is shown in the left panel of Fig. 8.The result in Eq. ( 45) is in good agreement with the estimate based on using the n = 0 screening mass and the dispersion relation, Eq. ( 44).
Using the same procedure for the continuum extrapolation of the non-static screening mass as in the case of the static, we obtained the following result in the first Matsubara sector: E (ns) 0 /T = 9.57 (19) stat (11) sys . ( The continuum extrapolation is shown in the right panel of Fig. 8.
As a crosscheck, we also investigated the possibility of determining the non-static screening mass using the ratio of the non-static and static screening correlators.By making use of the following approximate formula, we fitted the ratio of the correlators using the ansatz on the right-hand side of this equation.
After performing the averaging over the obtained differences (E 0 )/T , using the window function with the same parameters as in the case of the static screening mass analysis, the resulting continuum extrapolation is shown in Fig. 9, right.The outcome of the averaging using the window function centered at 1.05 is shown for the three discretizations on the finest ensemble in Fig. 9, left panel.The continuum estimate for the gap between the non-static and static screening masses at n = 1: Adding this to the value of the static screening mass determined directly from fitting the data, Eq. ( 45), we obtain 8.81(10) stat (2) sys + 0.60( 12) stat (2) sys = 9.41 (16) tot .Adding it to the static screening mass estimate using the continuum dispersion relation, Eq. ( 44), we get 8.72(2) stat (1) sys + 0.60( 12) stat (2) sys = 9.32 (12) tot .Both results are in good agreement with the result of Eq. (46).

IV. RESULTS CONCERNING THE n = 2 MATSUBARA SECTOR
In this section, we summarize the determination of H E (ω 2 )/T 2 , i.e. in the second Matsubara sector, and of the difference H E (ω 2 ) − H E (ω 1 ) using the more general subtractions given in Eqs.(19) and (20).More detailed discussion and further results applying these subtractions can be found in Appendices B and C.

A. Continuum extrapolation of various estimators for H
As was discussed in Sec.II B, when determining H E (ω n ), one is not limited to subtract the static screening correlator with the same momentum as the non-static screening correlator, but other choices are also possible.We referred to the momentum of the static screening correlator by adding a subscript p and wrote H E,p (ω n ) for the estimator in Eq. ( 18).This subscript was not used in other sections of the paper, since we have subtracted the static screening correlator with p = ω 1 when determining H E in the first Matsubara sector in Secs.III A-III C. A more general subtraction given in Eq. ( 19), based on a general linear combination of the static screening correlators with coefficients given by α, is also possible.The results obtained in this way have the additional index α, which we also omitted in previous sections but reinstate in the following discussion.
We note that on a finite lattice, the values of H E,p,α (ω n ) with different p i and α i values could differ from each other, but the result for H E (ω n ) in the continuum limit estimated using H E,p,α (ω n ) with different p i and α i values should be the same.One can therefore explore various choices to reduce the lattice artefacts of the continuum extrapolation.
The exact form of the integrand is in Eq. ( 50).
The integrands using the general subtractions of Eq. ( 19) show very different behavior compared to the integrand obtained by using the standard subtraction with p 1 = ω 2 and α 1 = 1 (α i>1 = 0).When modelling the tail of the integrands, a slightly extended version of the modelling procedure of Sec.III B was used due to the different noise level of the screening correlators in the different sectors.We discuss these procedures in more detail in App. C. We restrict our examination to choices which are a linear combination of the static correlator with zero and non-zero momentum and with the following relation between the coefficients defined in Eq. ( 19) In App.C, we explore a wider range of α values, but here we provide the results for only three choices of parameters: (i) the standard subtraction with α = 1 and p 1 = ω 2 (ii) a subtraction with α = 2 and p 1 = ω 1 which has the largest cancellation when Taylorexpanding in x 3 among the various terms in Eq. ( 50) at short distances (iii) a subtraction with α = 3.5 and p 1 = ω 1 which has the smallest observed slope in the continuum extrapolation among the investigated values of α (see Tab. II of App.C) We note that Eq.50 uses the continuum notation.First, we discuss the simplest of these choices, the standard subtraction (i).For this case, the modelling of the static screening correlator at ω 2 is more challenging than it was in the first Matsubara sector, because of the worse signal-to-noise ratio above xT ∼ 0.8.The continuum extrapolation has a moderate but significant slope (left panel of Fig. 10), and the obtained continuum limit results scatter in a wide range as the histogram in the right panel of Fig. 10 shows.This is primarily due to the uncertainty in the modelling of  50) with α = 2, p 1 = ω 1 as well as α = 3.5, p 1 = ω 1 plotted with green and blue colors, respectively.The latter curve starts at around −40, which is not included on the plot.the static screening correlator at p 1 = ω 2 .The obtained central value -although having a large systematic error -is in slight contradiction with the physical expectation, namely that , which is a consequence of the positivity of the transverse channel spectral function, σ(ω).This ordering of H E s in various Matsubara sectors can be deduced from Eq. ( 12).The subtractions (ii) and (iii) involving the static screening correlator at momentum p = ω 1 and at vanishing momentum is advantageous over the one at ω 2 primarily due to three reasons.First, the data has better precision because it does not involve the static screening correlator at spatial momentum p = ω 2 , therefore one can start the modelling from a later point.The details of this modified modelling prescription are discussed in App. C. Second, the plateau region is more clearly pronounced than in the case of the static screening correlator in the second Matsubara sector, which results in smaller systematic errors.Furthermore, due to the freedom of choosing α, one can construct integrands which behave better, which are less weighted towards the noisy long-distance tail or have smaller cutoff effects.
As mentioned earlier, the integrand obtained using more general subtractions can be quite different from what we obtain using the standard subtraction.In the continuum, both G T st (ω n , x 3 ) and G T ns (ω n , x 3 ) are expected to have the leading singular behaviour with the same prefactor 4 .This explains in particular the observed smoothness of the integrand (50) for our standard choice.For subtraction (ii), the integrand would be proportional to |x 3 | at small |x 3 | in the vacuum; therefore, the thermal integrand remains finite in the continuum at x 3 = 0 (the difference of the static thermal and the vacuum correlator has been investigated in [49]).For all other values of α, the integrand with p 1 = ω 1 has a singular behaviour around x 3 = 0, even though the integral has a well-defined continuum limit.Fig. 11 shows this for the second Matsubara sector, where one can see that −h α=1,p 1 =ω 2 (ω 2 , x), -i.e. the integrand using the standard subtraction -starts with a finite, positive value at x = 0 and after having a modest peak it decays to zero, being quite noisy above xT ∼ 0.8.The integrands using alternative subtractions also go to zero at large distances as earlier but the decay is faster.At p 1 = ω 1 , the choice α = 2, results in a smooth integrand, while the one with α = 3.5, has a more singular behaviour in the vicinity of x 3 = 0: the latter integrand starts at x/a = 0 with a value of about −40 for the finest ensemble and the LL discretisation, then it changes sign at x/a = 1 and has a peak at x/a = 2, after which it decays to zero.Thus, there is a large cancellation among the point at x 3 ≲ 1/(4T ).We use the trapezoid integration rule also at short distances, as we have done earlier.
The resulting continuum extrapolations using choices (ii) and (iii) are shown in Fig. 12, left and right panel, respectively.As one can observe in Fig. 12, the cutoff effects can be markedly different using different subtractions.For instance, the subtraction (ii) leads to a huge cutoff effect (left panel of Fig. 12), the results at the coarsest ensemble even have a different sign than the continuum estimate.On the contrary, the subtraction (iii) has a very flat continuum extrapolation (right panel of Fig. 12).The slopes of the continuum extrapolations are also listed in Table II of App.C for these subtraction types.
The continuum estimates we obtained by applying the different subtractions of the type given in Eq. ( 19) are the following: For further details on the modelling and the other parameters that were used to obtain these continuum estimates, we refer to App. C. In addition, in App.C, we investigate a broader set of α values, from which we conclude that the continuum result obtained by having a flat continuum extrapolation with α = 3.5 is a result that is in good agreement with more or less all the other results (c.f.Fig. 20 and also Fig. 21).Therefore, we choose this value, H E (ω 2 )/T 2 = −0.90(10) stat (4) sys (55) as our final continuum estimate in the second Matsubara sector.

B. Direct evaluation of the difference
Using Eq. 20, we can directly evaluate the difference of H E s obtained in different Matsubara sectors, or more generally the difference H E (ω r ) − εH E (ω n ).The results obtained in this way can be compared to the results obtained by calculating H E (ω r ) and H E (ω n ) separately and thereby serve as a useful crosscheck of those results.Additionally, the difference is an interesting quantity for its own sake.For instance, the choice ε = 1 probes an integrand, which is non-singular and is only sensitive to photons at nonzero frequencies (see Fig. 1).
We calculated the linear combination H E (ω 2 ) − εH E (ω 1 ) in several ways and summarized the results in Table III of App. C.Among the results listed in Table III, we highlight here the ones which have the flattest continuum extrapolations.These are the subtractions that have ϵ = 1, α 1 = 11.25, and ϵ = 14, In both cases, the momentum of the subtracted non-static screening correlator is ω 1 and the momenta of the static screening correlators are p 1 = ω 1 and p 2 = 0, for the correlators multiplied by α 1 and α 2 , respectively.Thus, the integrand in the continuum has the exact form of The results for ε = 1 correspond to forming the difference of the integrand for H E (ω 2 ) that has a flattest observed continuum extrapolation (choice (iii) of the previous section, see Eq. 50) and the integrand for H E (ω 1 ) having a standard subtraction term (Eq.18 with n = 1 and p = ω 1 ).Therefore, it is nice to observe that the difference obtained in the continuum limit 1 is in complete agreement with the difference formed by using the continuum estimates for H E (ω 2 ) and H E (ω 1 ), separately.
For the other set of parameters with ϵ = 14 for which we also observe a mild continuum extrapolation, we obtain p=(ω 1 ,0), α=(0,−13) = 8.44 (11) stat (5) sys . ( In order to put this into context, we calculated H E (ω 2 )/T 2 from this result by adding 14 times the continuum estimate of H E (ω 1 )/T 2 in a correlated way.We obtained H E (ω 2 )/T 2 = −0.94(12) stat (5) sys , which is also in good agreement with our final estimate.The integrand for determining H E (ω 1 ) in the free theory (dashed) and using lattice QCD data on our finest ensemble X7, completed with a model based on single-state fits for long distances using Eqs.( 30),( 34),(38) (solid line).Only the central values are shown.

V. COMPARISONS
In this section, we compare our findings to results from analytic approaches.We start with the integrand for H E (ω 1 ), which is shown in Fig. 13.The calculation of the integrand of the free continuum result requires special care at very short distances, below xT ∼ 0.1.This is discussed in App.D. By looking at Fig. 13, we can observe that the continuum free theory result has very different characteristics compared to the lattice QCD result.Although both start with a finite value at x = 0, the free theory result decays faster and even goes to slightly negative values above xT ≈ 1.
Secondly, we summarized the comparison of the non-static screening masses determined on the lattice and the weak-coupling results in Fig. 14.We refer to App.F as well as to Ref. [37] for more detailed information about the results.As Fig. 14 shows, the lattice result for the non-static screening mass in the n = 1 sector lies between the NLO and the EQCD result.We conclude that the weak-coupling prediction is fairly successful, at the 4% level, provided at least the next-to-leading order interquark potential is used.The weak-coupling results in Fig. 14 are based on the value α s = 0.25, which will be our default value of the gauge coupling in the following.
Calculating the imaginary part of the retarded correlator at lightlike kinematics in the free massless theory with our current normalization gives H E (ω n )/T 2 = −N c /6, independently of n [31].This corresponds to a free spectral function Therefore H E (ω 1 )/χ free s = −0.500.This result can also be reproduced to good precision by integrating the free theory integrand shown in Fig. 13, if one uses a suitable representation (cf.appendix D).
In strongly coupled N = 4 super Yang-Mills (SYM) theory using the AdS/CFT correspondence, one finds H E (ω 1 )/χ free s = −0.336[31], using (χ s ) free = N 2 c T 2 /4 (see [51], appendix A).This value is in fact lower than the free-theory result.Compared in this way, the lattice result we obtained, −0.670( 6) stat (1) sys , is largest.Using a different normalization, e.g.dividing by the static susceptibility, χ s , of the relevant interacting theories we arrive at the following predictions: in the free theory, obviously the result does not change, while in N = 4 super Yang-Mills theory, we get [H E /χ s ] (SYM) = −0.6715[31], where χ s = N 2 c T 2 /8.In Ref. [23], we determined the static susceptibility at this temperature to be χ (lat) s /T 2 = 0.882 (11) stat (19) sys .(62) Using this value, our lattice result is [H E /χ s ] (lat) ≃ −0.76.Thus, now normalizing by the interacting χ s , the lattice result is still the largest in magnitude.These results are illustrated in Fig. 15.
In the second Matsubara sector, we have −1.115 for H E (ω 2 )/T 2 in strongly coupled N = 4 SYM [31], thus the difference (H E (ω 2 ) − H E (ω 1 ))/χ s is −0.444, a value larger than our lattice result by 1.7 standard deviations.Since H E (ω n ) is constant in the free theory, the difference vanishes there.Thus the ratio (H E (ω 2 ) − H E (ω 1 ))/H E (ω 1 ) provides good sensitivity to the shape of σ(ω), being 0.66 in the strongly-coupled SYM case and parametrically small in the case of a spectral function very peaked around ω = 0.The lattice result, (H E (ω 2 ) − H E (ω 1 ))/H E (ω 1 ) = 0.34 ± 0.19, lies between these two extremes.

A. Computing H E (ω n ) dispersively using the complete leading-oder σ(ω)
Using the complete leading-order result of Arnold, Moore and Yaffe (AMY) for σ(ω) [12,13], the difference [H E (ω 2 ) − H E (ω 1 )] can be evaluated straightforwardly, whereas the quantity H E in individual Matsubara sectors can only be estimated after handling the singular behavior of σ AMY (ω) at small frequencies.Indeed, by integrating the spectral function of Ref. [12,13] (with α s = 0.25) multiplied by the kernel ω 2 n /(πω(ω 2 + ω 2 1 )) in the range 0.2 < ω/T < 50 where the provided parametrisation is a good approximation, we obtain which is already larger in absolute value than the result obtained on the lattice.On the other hand, the prediction for H E (ω 2 ) − H E (ω 1 ) using the AMY spectral function is much less sensitive to the small-ω behaviour.We obtain, again with α s = 0.25, The comparison can be seen in Fig. 15, right panel.The leading-order prediction for this difference is well compatible with our lattice QCD result.We also note that extending the integral to ω = ∞ assuming σ(ω) ∝ √ ω [12] changes the value of Eq. (64) to −0.28.In order to exploit our precise lattice result for H E (ω 1 ), we need to inspect more precisely the region of validity of the weak-coupling spectral function.The leading-order calculation [12,13] assumes the photon wavelength to be short compared to the mean free path for large-angle scattering.At the smallest frequencies at which the calculation is still valid, σ(ω) ∝ 1/ √ ω.Since we expect σ(ω)/(2χ s ω) to tend to a finite value at ω → 0, namely to the diffusion coefficient D (see Eq. ( 14)), it is clear that a qualitative change in the functional form must take place below the frequency at which the AMY calculation breaks down.In addition, the next-to-leading order correction [52], suppressed only by one power of the strong coupling g, turns out to be quite modest (about ±10%) for ω ≳ 2πT , but becomes larger (about ±30%) at ω ≲ πT .
Thus a qualitative modification of the AMY spectral function at small frequencies is necessary for a sensible dispersive evaluation of H E (ω 1 ).In the following, we assume that σ(ω) is given by the leading-order expression [12] for ω > ω m , introduce a Lorentzian ansatz for ω < ω m , σ(ω) Beyond a value ω = ω m marked by a short vertical line, the curves correspond to the complete leading-order result in the parametrization provided by [12], for two values of α s .For ω < ω m , the curves correspond to the Lorentzian (65) with parameters given in Eqs.(66-67), respectively.In the α s = 0.25 case, the dashed curve shows how the parametrization of [12] extends towards smaller frequencies.The intercept at ω = 0 yields the isospin diffusion coefficient, T • D. The leading-order weak-coupling result for D is known from Ref. [53]: with m D /T = 2.05, which results from setting α s = 0.25 in the leading order expression of the Debye mass, we read off T • D = 2.3 from Fig. 1 of [53].It turns out that this value of the diffusion coefficient can be reproduced by choosing the matching point at ω m = 0.50T , using α s = 0.25 as well in the weak-coupling spectral function σ(ω) [12].In that case, however, the dispersive integral yields H E (ω 1 )/T 2 = −0.99, in stark disagreement with our lattice result (41).Thus this particular weak-coupling scenario is ruled out at T = 250 MeV by our lattice calculation.
Instead of assuming the weak-coupling result for the diffusion coefficient, we can attempt to estimate it by continuing the leading-order expression [12] of σ(ω) toward the soft-photon limit via the Lorentzian Eq. ( 65) so as to reproduce our lattice result for H E (ω 1 ).With α s = 0.25, this condition yields Increasing the targeted |H E (ω 1 )|/T 2 by one standard deviation (0.006) only results in a modest change in the estimated diffusion coefficient to T • D = 0.56.Thus the estimate of the diffusion coefficient obtained in this way is much lower than the weak-coupling result [53], but still more than three times larger than the AdS/CFT value of D = (2πT ) −1 .Repeating the procedure above with a larger value of the coupling, we obtain In summary, under the assumptions made on the spectral function, we obtain a range of diffusion coefficients given by Eqs.(66-67).The corresponding spectral functions are illustrated in Fig. 16.
It is worth recalling that in the non-interacting theory, the spectral function has the form ωδ(ω) (see Eq. ( 61)).In the limit of vanishing coupling, one thus expects the Lorentzian to turn into a delta function (η → 0 at fixed A/η), and the spectral function to vanish roughly proportionally to α s for ω ≫ η.We define the area under the peak around ω = 0 as follows, where Ω ≫ η is a UV-cutoff.On the one hand, we get with our Lorentzian ansatz S 0 = A/η.In the free theory, one obtains S 0 = 1/2.Thus it is remarkable that the ratios A/η in Eqs.
(66-67) are close to 1/2: the obtained values of the parameters (A, η) are plausible in this respect, since our analysis is based on an ansatz for the spectral function valid at weak coupling and S 0 has a weak-coupling expansion of which S 0 = 1/2 is the leading term.

VI. CONCLUSIONS AND OUTLOOK
The thermal photon emissivity of the quark-gluon plasma is determined to all orders in the strong coupling by the transverse channel spectral function of electromagnetic currentcurrent correlators evaluated at lightlike kinematics.This real-time observable is not accessible directly on a Euclidean lattice.In this work, we demonstrated that it is nonetheless possible to directly evaluate an observable in lattice QCD that is related to the aforementioned spectral function via a dispersion relation.Indeed, the computed observable -a Euclidean screening correlator at imaginary spatial momentum -has an integral representation in terms of a product of the spectral function multiplied by a Lorentzian kernel (see Eq. 13).However, the naive lattice estimator of this observable is afflicted by a large cutoff effect arising from the breaking of Lorentz invariance on the lattice.We successfully addressed this technical difficulty by subtracting screening correlators at different momenta, whose contribution vanishes in the continuum.We investigated various subtractions which lead to different scaling behaviors towards the continuum limit.
By analyzing the long-distance behavior of the screening correlators, we also determined the screening masses (see Eqs. (43,45,46)), which we used to model the position-space correlators at long distances and improve the precision of the imaginary-momentum correlators.In order to probe the virtuality dependence of the latter around the Q 2 = 0 point, we constructed a suitable lattice representation of its Q 2 -derivative and evaluated it in the first Matsubara sector.We used two-flavors of O(a)-improved Wilson fermions at a temperature of about 250 MeV with an in vacuo pion mass of 270 MeV and performed continuum extrapolations of all our observables.
Our final continuum result for H E (ω 1 )/T 2 , see Eq. ( 41), has a precision of around 1%, suitable for a comparison with other approaches.We confronted our results to estimates using the free theory, strongly-coupled N = 4 SYM and the full leading-order result of Arnold, Moore and Yaffe (AMY) [12,13].We found that our result for |H E (ω 1 )| is smaller than the prediction obtained by using the spectral function of Ref. [12,13], although our result for H E (ω 2 ) − H E (ω 1 ), given in Eq. ( 59), is comparable.This is illustrated in Fig. 15.Since the integration kernel is much more sensitive to the low-frequency behavior of the spectral function in the case of H E (ω n ) than in the case of the difference (c.f.Fig. 1), our result suggests that the weak-coupling result for the spectral function is overestimated in this low-frequency region.Assuming the AMY spectral function to hold above a certain frequency ω m , and using a Lorentzian transport peak below that frequency that matches on smoothly at ω m , we arrive at estimates of the isospin diffusion coefficient T • D in the range 0.35 to 0.54 (see Eqs. 66-67) by requiring that our lattice result for H E (ω 1 ) be reproduced.This range of values for T • D is in line with previous lattice estimates based on dispersion relations at fixed spatial momentum [20,25,38,[54][55][56][57], though the central values of most calculations with dynamical quarks lie below the value of 0.3 at temperatures around 250 MeV [17].
It is also interesting to compare our results with our two previous studies of the photon emissivity that employed the same gauge ensembles as the present calculation but were based on the dispersion relation at fixed spatial momentum [23,49].Addressing the inverse problem with physically motivated ansätze for the spectral functions, we concluded [23] (particularly for ω ≥ πT ) that the lattice results were consistent with the weak-coupling prediction, but could also accomodate a rate 2.5 times larger.In other words, most of the solutions for the spectral function describing the lattice data yielded a photon rate at least as large as the weak-coupling prediction.Given that the weak-coupling spectral function [12] results in a larger value for |H E (ω 1 )| than our lattice result, it seems most likely that this excess is due to an overestimated soft-photon emissivity in the weak-coupling calculation.
We have seen that computing H E (ω 2 ) by standard lattice methods is already a lot more challenging than for H E (ω 1 ) due to the large statistical errors on the position-space integrand.Thus it would be interesting to investigate noise-reduction methods similar to those used in calculations of the hadronic vacuum polarization.Improving on our determination of H E (ω 2 ) − H E (ω 1 ) would go a long way to ascertain that the emission rate of hard photons (ω ≳ πT ) follows the weak-coupling prediction at the temperatures reached in heavy-ion collisions.This appears to be an achievable goal in the near future.
A second difficulty that occurred is the increase in discretization errors for increasing ω n .Perhaps some improvements in the discretization are still possible here, but it is clear that very fine lattices are necessary in order to determine H E (ω 3 ) and beyond.While the moments H E (ω n ) provide valuable non-perturbative constraints on the spectral function, a numerically ill-posed inverse problem resurfaces if one has the ambition of determining σ(ω) itself from a (necessarily) finite collection of its moments.Instead, our immediate plan is to compute the first moment of σ(ω) across the phase crossover with physical quark masses in order to test the models of σ(ω) used in hydrodynamics-based calculations of photon spectra in heavy-ion collisions.trum Jülich allocated under NIC project HMZ21.For generating the configurations and performing measurements, we used the openQCD [45] as well as the QDP++ packages [58], respectively.

Appendix A: Elimination of outliers
In this Appendix, we describe our procedure to deal with certain measurements that deviate by several standard deviations from the mean value calculated using the data.We identify these exceptional measurements as outliers, see the left panel of Fig. 17 and we found, they occur more frequently at large Euclidean separations.These outliers increased the statistical error and also modified the mean to some extent, see Fig. 17, right panel.We eliminated them by using robust statistics [59].
In our procedure, we first prepared a distribution of results at each Euclidean distance, then removed the data points belonging to the lower and upper γ% of that distribution.Varying γ in the interval 0.5-4, we made cuts and found that the error estimation as well as the calculation of the mean is more stable this way.When we detected only less than 10 datapoints being outside five times the interquantile range from the mean, we applied only a trimming with γ = 0.5.For other distances we applied γ = 1 in our final analysis.
We show an example in the case of the conserved-conserved correlator on our finest ensemble, X7, in the right panel of Fig. 17.At short distances of the correlator, this approach did not influence the results, because outliers occured there only very rarely.At intermediate distances, i.e. around xT ∼ 0.7-1.3, the effect of this method was again not significant.At large distances, however, the errors reduced by a factor of around 2-6 when omitting the tails of the distributions.Although one has to be careful when discarding certain measurements, we believe that our procedure removing outliers sometimes more than 10 standard deviations off from the mean at large distances should not influence the validity of the extracted physical results.When subtracting the static screening correlator at p = ω 1 or at p = 0, we slightly modified the procedure of modelling the integrand.We divided the integration interval into three parts.At short distances, both the non-static and static screening correlators have good signal-to-noise ratios, therefore we applied no modelling there.At the intermediate interval, the non-static screening correlator is already modelled but the static screening correlator is not.In this interval, we calculate the integral using the trapezoidal rule applied for the integrand formed by evaluating the single-state fit used for modelling the non-static screening correlator at the lattice points and subtracting the lattice data for the static screening correlators.Finally, in the interval for large distances, we modelled the static screening correlators as well and applied Eqs.(38) and (39).We start applying modelling by single-state fits using the smooth step function of Eq. ( 35) from x w,ns T = 0.6, 0.7 in the case of the non-static screening correlator in the second Matsubara sector and from x w,st T = 1.0, 1.1, 1.2, 1.3 in the case of the static screening correlator at p = ω 1 or at p = 0.
After evaluating the integrals, we performed the continuum extrapolations in a similar manner as was discussed in Sec.III C. The continuum limit fit form is where A 0 is the estimate of the continuum limit from a particular fit, and the A δ 1 parameters characterize the approach to the continuum of the H E values calculated using the different discretized correlators.δ stands for LL, LC or CC.
First, we present results that were obtained using the simple alternative subtraction of Eq. ( 18) but with p = ω 1 or p = 0.Although the lattice artefacts are huge in the case of H E,p=0 (ω 2 )/T 2 and H E,p=ω 1 (ω 2 )/T 2 , the fit qualities of the continuum extrapolations using Eq.(C1), i.e. a linear fit ansätze in a 2 turned out to be acceptable.For instance, only around ∼ 2% of all the p-values is smaller than 0.05.The estimates in the continuum using these subtractions are: H E,p=0 (ω 2 ; a)/T 2 = −0.74(8) stat (5) sys , (C2) lim H E,p=ω 1 (ω 2 ; a)/T 2 = −0.71(12) stat (4) sys . (C3) Depending on the actual modelling interval, the results could be slightly different, but stay consistent within errors.These subtractions bring only a modest improvement: the modelling is more precise, but the slope of the continuum extrapolations are huge (see Fig. 19).
Applying the more general subtractions of Eq. ( 19) and Eq. ( 20), we determined H E (ω 2 )/T 2 using a broad set of parameters.We discussed a subset of these results in Secs.IV A and IV B. The results employing a subtraction with vanishing ε are labelled as fitI and are listed in Table II.For those subtractions we used a single parameter, α, with which the integrand can be described as in Eq. (50).Besides the parameters that characterize the start of the modelling of the non-static as well as the static screening correlators using single state fits (T w,ns and T x w,st , respectively), we also list the fit parameters (A LL 1 , A LC 1 , A CC 1 ) that characterize the slopes of the continuum extrapolations in Table II.Using these values, one can read off that the subtraction of the form of Eq.50 gives the flattest continuum extrapolation with the choice of α = 3.5.With this value, the slope parameters of the continuum extrapolation are consistent with zero for the local-conserved and conserved-conserved discretizations.We summarize the various results for H E (ω 2 )/T 2 using these type of alternative subtractions in Fig. 20.
Turning to the subtractions that involve the non-static correlators at ω 2 and ω 1 as well, we applied a similar 'two-interval' modelling procedure as we discussed above, but used a larger T x w,st value for the non-static screening correlator in the first Matsubara sector as well.The results employing this type of subtraction with non-vanishing ε are labelled as fitII and are listed in Table III.In Table III, all results correspond to the choice (T x w,ns , T x w,st ) = (0.7, 1.1); as we did earlier (see e.g. in Table II), we performed a scan varying these parameters, but found only small changes.
It is useful to recall that the subtractions using Eq. ( 20) (see also Eq. ( 58)) enable a direct calculation of the difference of H E s in different Matsubara sectors.We discussed two of these type of results in Sec.IV B. In Table III, more of these types of results are reviewed.Since the continuum extrapolated quantity is [H E (ω 2 ) − εH E (ω 1 )] p,α in this case, we also included a column denoted by A 0 , that explicitly contains this estimate.By adding ε • H E (ω 1 )/T 2 to this value, one can obtain an estimate for H E (ω 2 )/T 2 .We did this in a correlated way, and then estimated the statistical error on H E (ω 2 )/T 2 and added the systematic errors in quadrature.The results of Table III (35) Table II: Results for H E (ω 2 )/T 2 using alternative subtractions of the form given in Eq. (19).The values of the momentum in the subtraction given in Eq. ( 19) are p 1 = ω 1 and p 2 = 0. α is a single parameter in Eq. ( 19), that fixes α 1 and α 2 : propagator.Consider the function

Now we can write
(where we have included the color factor N c = 3 explicitly) and The integrand for H E (ω n , Q 2 = 0), defined to be negative-definite, is given by the difference of the two preceding correlators: Using the representation (D3) of the function H, the 'no-winding' term n s = 0 cancels in the difference G T ns (ω n , x 3 ) − G T st (ω n , x 3 ); this cancellation corresponds to the fact that the integrand vanishes in the vacuum.The integrand can then be evaluated efficiently at small x 3 , even directly at x 3 = 0.At x 3 ≳ 1/(2πT ) on the other hand, the expressions (D5) and (D8) should be used to evaluate the correlator as a rapidly converging sum.

E
The free Wilson quark propagator is diagonal in color space and can be written in the time-momentum representation as ⟨ψ(x) ψ(y)⟩ with the natural convention sgn(0) = 0 for the sign function and the standard notation The constants appearing in the numerator are given by The single-quark energy pole is given by Further constants are The denominator reads After these preliminaries, we are ready to compute the lattice screening correlators at the one loop level.For that purpose, we use x 3 as the 'time' direction, so that the quark correlator falls off like exp(−ω ⃗ p |x 3 |).Here we give the result for the local-conserved discretization; the expression of the local and the conserved currents are given in Eqs. ( 25) and ( 26) respectively.Let ⃗ k = (k 0 = ω n , k 1 , k 2 ) be the external momentum, ω n = 2πT n, n = 0, 1, . . .N t − 1. Define With ⃗ p = (p 0 , p 1 , p 2 ) = ((2ν + 1)πT, p 1 , p 2 ), ν = 0, . . ., N t − 1, the free-quark prediction is where N c = 3 is the color factor and We are interested in the massless theory, in which case no renormalization factor or additive O(a)-improvement is needed at leading order in perturbation theory, since in continuum perturbation theory the vector-tensor correlator vanishes in the chiral limit.Thus we expect G CL 11 ( ⃗ k; x 3 ) to be O(a)-improved at a fixed x 3 ̸ = 0.In the free massless theory, the integral over x 3 leading to H E does not necessarily converge at long distances, due to the singlequark lattice dispersion relation being modified by O(a 2 ) from its continuum counterpart.Therefore, in order to judge the size of cutoff effects, we consider in each Matsubara sector the (discretized) truncated integrals with The expressions above correspond to choosing the trapezoidal rule for the corresponding integrals.In addition, we consider the estimator The static contributions appearing in H II E and H III E do not vanish as they would if the integral extended to x 3 = ∞.Therefore one should not expect H I E , H II E and H III E to have the same continuum limit.However, in the interacting theory, we expect the bulk of the discretization errors to come from the region 0 ≤ x 3 ≤ β.Each of the three quantities leads to a separate estimator for H E by extending the integral to x 3 = ∞.We therefore investigate the cutoff effects on H I E , H II E and H III E in order to assess the relative merits of the corresponding estimators for H E .The remarks around Eq. ( 51) concerning the behaviour of the integrands in the vicinity of x 3 = 0 for the various estimators apply in particular to the free case investigated in this appendix.
Table IV  Table V provides a similar comparison in the Matsubara sector ω 2 .As one might expect, cutoff effects are overall larger in this sector.The quantity H II E approaches its continuum limit much faster than H I E , though clearly, if data were only available up to N t = 24, an extrapolation would be needed to reach a precision of a few percent on the continuum result.In the lattice QCD simulations, however, the static correlator with spatial momentum equal to ω 2 tends to be noisy.This motivates us to consider alternative subtraction schemes such as H III E (α, ω 2 ).The choice of α = 2, while providing a smooth x 3 integrand, does not improve the approach to the continuum as compared to H I E .On the other hand, a somewhat larger value of α (for instance 3.5) does lead to reduced cutoff effects.It is worth pointing out that, in the continuum, the integrals for H I,II,III E all converge to −0.5 when extended to x 3 = ∞.The quantity H III E (3.5, ω 2 ) is already quite close to that value, indicating that the integrand in this case must be very suppressed for x 3 > β, which is a further desirable feature.

Appendix F: Weak-coupling theory results for the screening masses
In this appendix, we overview our results for the screening masses obtained using weakcoupling theory.
As it has been derived in Ref. [37], in order to calculate the spectrum of non-static screening states in weak-coupling theory, one has to solve an inhomogeneous Schrödinger equation.Interestingly, apart from different normalizations, this has the same general form as the one in the LPM resummation of the photon production rate.By solving the radial part of the homogeneous Schrödinger equation numerically, we can determine the non-static screening energies by plugging in the obtained eigenvalues, Ê(l=1) , into Here, The value of the gauge coupling of the dimensionally reduced effective theory at two loops is determined to be g  Table VIII: The non-static screening energies using the EQCD potential determined at T = 400 MeV [62].
Appendix G: One-derivative operators having an improved overlap onto the ground state in the transverse non-static sector In order to consolidate our extraction of the screening masses from the transverse-channel non-static screening correlators, Eq. ( 27), we look for alternative operators that may overlap better with the low-energy states.We therefore consider non-static two-point functions of quark bilinears in the free theory at infinite spatial volume: with s 0 = sgn(x 0 − w 0 ).Time is in the direction "0".
There are four of such one-derivative operators, which give a non-vanishing Dirac trace when correlated with the current ψγ 1 ψ.However, when considering the C-parity of these operators, one finds that the ones in Eq. (G5) are C-even.Recalling that the conserved vector current is C-odd, we find that correlation function of the operators in Eq. (G5) with ψγ 1 ψ vanishes.Therefore only the operators of Eq. (G4) need to be considered.The Dirac traces T (Γ, Γ) of these operators with plus (minus) sign for the operator with Γ = Γ = γ 0 γ 5 (Γ = Γ = γ 3 γ 5 ).Thus, we find that in the free theory both operators of Eq. (G4) have an unsuppressed coupling to low-lying states.
We have investigated the two-point function of the second operator of Eq. (G4), now in the interacting theory on our coarsest ensemble called F7.This 12 × 48 3 ensemble at the same temperature of T = 254 MeV has not been included in the main analysis.Its relevant parameters are given in Table I of Ref. [23].
We introduce the following notation for the lattice operator, where the forward-backward differential  Using the above definitions, we then calculated the two-point correlation function of the O 35D (x) operator.We performed this evaluation using also a smeared operator, which reduced the noise significantly.We also note that, similarly to the transverse screening correlators of the vector current, we performed averages over different spatial decay directions exploiting the discrete lattice rotation symmetry.
We extracted the local effective mass from the two-point function and we also performed single-state fits using a similar fit ansatz as in Eq. (32).Using the fit results for the mass, we constructed the window-smeared mass averaging fit results obtained at subsequent equallength fit ranges.The window function was centered at xT = 0.6 and had a width of 0.30 in this variable.These correspond to x win,center /a ≈ 7 and a width of 3.6 in lattice units.We note that for instance the fit results using 9 datapoints starting with x/a = 4, 5, 6, 7, 8, 9, 10 have p-values 0.02, 0.27, 0.15, 0.18, 0.45, 0.68, 0.59, respectively, i.e. the chosen window range covers results with acceptable p-values.The window-averaged mass is shown by the horizontal blue bland in Fig. 22.Compared to this we also show in Fig. 22 the singlestate fit results for the two-point function involving the smeared or unsmeared V 1 operator.The window-smeared mass obtained using the unsmeared V 1 is shown with the short red horizontal band in Fig. 22 and is in good agreement with the one obtained using O 35D (x).
In summary, while the smeared O 35D and V 1 correlators exhibit lower effective masses in the region x 3 T ≲ 0.75, their variance is also larger.This observation should be taken into account in future dedicated spectroscopic calculations.

Figure 2 :
Figure2: The integrand (Eq.(30)) for the calculation of H E (ω n ) on our finest ensemble X7, using the conserved-conserved discretization of the current-current correlator.

Figure 3 :
Figure 3: Left: fitted masses, with circles and squares corresponding to G T ns (ω 1 , x 3) and G T st (ω 1 , x 3 ) respectively, using a fit range of ten lattice spacings.The bands represent the 'cosh' effective masses.Right: the tail of the integrand needed for the calculation of H E at the first Matsubara frequency.The red points are the actual datapoints on our second coarsest lattice and the band shows the result of the modelling.

Figure 4 :
Figure4: Left: A representative continuum limit of −H E (ω 1 )/T 2 obtained from Eqs.(29)(30).Right: AIC-weighted histogram of the continuum extrapolated results for −H E (ω 1 )/T 2 obtained with x w T = 1.1.The AIC-weighted histogram of the long-distance contribution to H E (ω 1 ), shifted with the continuum result for the short-distance contribution, is also shown for comparison (sd+ld).The dark-grey band shows the systematic error, and the lighter gray band represents the total error obtained from the statistical and systematic errors added in quadrature.

Figure 5 :
Figure 5: Comparison of continuum results for −H E (ω 1 )/T 2 using different starting points for the modelling.

Figure 7 :
Figure 7: Continuum extrapolation of the static screening mass obtained using the screening correlator at vanishing momentum.

Figure 9 :
Figure 9: Left: Fit results for the gap between the ground-state non-static and static screening masses, respectively associated with the correlators G T ns (ω 1 , x 3 ) G T st (ω 1 , x 3 ), and the windowsmeared results shown by the bands.Right: Continuum extrapolation of the difference of the ground-state non-static and static screening masses.

Figure 11 :
Figure11: Comparison of the integrands obtained by applying the standard subtraction (i) (red), and more general subtractions (ii) and (iii) of the form given in Eq. (50) with α = 2, p 1 = ω 1 as well as α = 3.5, p 1 = ω 1 plotted with green and blue colors, respectively.The latter curve starts at around −40, which is not included on the plot.

Figure 13 :
Figure13: The integrand for determining H E (ω 1 ) in the free theory (dashed) and using lattice QCD data on our finest ensemble X7, completed with a model based on single-state fits for long distances using Eqs.(30),(34),(38) (solid line).Only the central values are shown.

Figure 14 :
Figure14: Ground-state non-static screening masses in the n = 1 Matsubara sector determined on the lattice and by using weak-coupling theory at leading order (LO), with an interquark potential at next-to-leading order (NLO), and by using the potential at T = 400 MeV of the dimensionally reduced effective theory for QCD, 'electrostatic QCD' (EQCD).

31 Figure 16 :
Figure16: Two weak-coupling based spectral functions reproducing the value(41) for H E (ω 1 ).Beyond a value ω = ω m marked by a short vertical line, the curves correspond to the complete leading-order result in the parametrization provided by[12], for two values of α s .For ω < ω m , the curves correspond to the Lorentzian (65) with parameters given in Eqs.(66-67), respectively.In the α s = 0.25 case, the dashed curve shows how the parametrization of[12] extends towards smaller frequencies.The intercept at ω = 0 yields the isospin diffusion coefficient, T • D.

Figure 17 :
Figure 17: Left: Measurement history of the non-static screening correlator at xT = 0.833.The outliers are shown up as spikes in the data.Right: Truncation stability for the non-static screening correlator.The data points corresponding to the trimmed data have been shifted slightly to the right to improve visibility.

Figure 22 :
Figure22: Fitted and local effective (cosh) masses with symbols and bands, respectively, extracted from the two-point correlation functions of the V 1 , smeared V 1 and smeared O 35D operators with red circles, green squares and blue triangles, respectively.The short, horizontal red and blue bands show the value and error of the window-averaged masses.
Figure15: Left: Comparison of results for H E (ω 1 ) obtained on the lattice, in leading order of the weak-coupling expansion (with α s = 0.25 (α s = 0.31) for the left (right) point), in the free theory and in strongly coupled N = 4 SYM, normalized by the static susceptibility of the corresponding non-interacting theory (left inset), or by the interacting static susceptibility (right inset).We interpret the two normalizations as being identical for the "free" data points, while to normalize the NLO results we used alternatively 1 and the lattice result 0.88 for χ s /T 2 .Right: Comparison for H E (ω 2 ) − H E (ω 1 )with the same choices of normalization.
are summarized for better overview in Fig.21. Figure 20: Comparison of continuum extrapolations for H E (ω 2 )/T 2 using various subtractions.The parameters of the various fits are listed inTable II.The final estimate for H E (ω 2 )/T 2 is shown with a blue vertical band.The final estimate for H E (ω 1 )/T 2 is also shown for comparison with a green vertical band.fit α T x w,ns T x w,st −H E (ω 2 )/T 2

Table III :
(19)re21: Comparison of continuum extrapolations for H E (ω 2 )/T 2 using various subtractions involving also the non-static correlator at ω 1 .The parameters of the various fits are listed in TableIII.The final estimate for H E (ω 2 )/T 2 is shown with a blue vertical band.The final estimate for H E (ω 1 )/T 2 is also shown for comparison with a green vertical band.Results for H E (ω 2 )/T 2 using alternative subtractions with ε ̸ = 0.The formula for the integrand is given in Eq.(19).The continuum limit fit ansatz parameters (see Eq. C1) are also listed in the table.

Table IV :
Approach to the continuum of different integrals for the first non-zero Matsubara sector ω 1 .

Table V :
II E .It is likely related to the fact that the corresponding estimator for H E (ω 1 ) is free of cutoff effects in the vacuum, and therefore any cutoff effect must depend mainly on the parameter aT = 1/N t .N t H I E (ω 2 )/T 2 H II E (ω 2 )/T 2 H III E (3.5, ω 2 )/T 2 H III E (2.0, ω 2 )/T 2Approach to the continuum of different integrals for the second non-zero Matsubara sector ω 2 .
compares the approach to the continuum for the two quantities H I E and H II E in the Matsubara sector ω 1 .Clearly, the approach is much faster for the quantity H

Table VI :
The numerically determined eigenvalues as well as the corresponding energies for the LO, NLO and EQCD cases are listed in the Tables VI, VII and VIII, respectively.The non-static screening energies at LO.

Table VII :
The non-static screening energies using the NLO potential.