Heavy Quark Diffusion from 2+1 Flavor Lattice QCD with 320 MeV Pion Mass

We present the first calculations of the heavy flavor diffusion coefficient using lattice QCD with light dynamical quarks. For temperatures $195\,\mathrm{MeV}<T<352\,\mathrm{MeV}$, the heavy quark spatial diffusion coefficient is found to be significantly smaller than previous quenched lattice QCD and recent phenomenological estimates. The result implies very fast hydrodynamization of heavy quarks in the quark-gluon plasma created during ultrarelativistic heavy-ion collision experiments.

Introduction.-Heavy charm and bottom quarks are produced only during the earliest stages of ultrarelativistic heavy-ion collisions. They participate in the entire evolution of the quark-gluon plasma (QGP), and emerge as open heavy-flavor hadrons or quarkonia. These heavyflavor hadrons provide valuable insights into the QGP. Experiments at the relativistic heavy ion collider (RHIC) and large hadron collider (LHC) show that the yields of heavy flavor hadrons at low transverse momentum (p T ) are asymmetric along the azimuthal angle, meaning that the elliptic flow parameter (v 2 ) is large [1][2][3]. These observations indicate that heavy quarks participate in the hydrodynamic expansion of the QGP. Understanding the origin of the hydrodynamic behavior of heavy quarks is key to understanding the near-perfect fluidity of the QGP.
The motion of a heavy quark with mass M T , immersed in a QGP at temperature T , can be effectively described by Langevin dynamics [4]. The heavy quark momentum diffusion coefficient, κ, quantifies the momentum transfer to the heavy quark from the QGP background through random momentum kicks which are uncorrelated in time. Specifically, 3κ is the mean squared momentum transfer per unit time. For sufficiently low p T , hydrodynamization of heavy quarks in the QGP can be characterized by the diffusion constant in space D s = 2T 2 /κ [4], where 6D s is the mean squared distance traversed per unit time by the heavy quark inside the QGP. The available perturbative QCD results for κ [5] provide reliable estimates only for asymptotically large T . For practical phenomenological studies of heavy-ion collisions lattice QCD results for κ are needed. Until now, lattice QCD based determinations of κ have been limited only to pure gauge theory, therefore neglecting dynamical fermions entirely. Here, we report the first continuum-extrapolated lattice QCD calculations of κ in 2+1-flavor QCD with a physical strange quark mass and degenerate up and down quark masses corresponding to a pion mass m π 320 MeV.
Theoretical framework.-The momentum diffusion coefficient κ can be obtained from the zero-frequency limit of the spectral function corresponding to the conserved current-current correlation function of heavy quarks. Relying on M T and M Λ QCD , the QCD Lagrangian can be expanded in powers of 1/M . By integrating out the heavy quark fields one arrives at the heavy quark effective theory. In this effective theory the heavy-quark current-current correlator at leading order in 1/M becomes equivalent to the correlation function of the colorelectric field [6,7]: ReTr [U (1/T, τ )E i (τ, 0)U (τ, 0)E i (0, 0)] 3 ReTrU (1/T, 0) .
Lattice QCD calculations.-We performed calculations in 2+1 flavor QCD with a physical strange quark mass, m s , and degenerate up, down quark masses m l = m s /5 using the highly improved staggered Quark (HISQ) action [9] and tree-level improved Lüscher-Weisz gauge action [10,11]. In the continuum limit our choice of m l corresponds to m π 320 MeV. The lattice spacing a and the quark masses are fixed as in Refs. [12,13]. We carried out calculations on 96 3 × N τ lattices with 1/a = 7.036 GeV and N τ = 20, 24, 28, 32 and 36, that correspond to temperatures T = 352, 293, 251, 220 and 195 MeV, respectively. To control discretization effects we also performed calculations on 64 3 × N τ lattices (N τ = 20, 22 and 24) at different lattice spacings, chosen such that the above temperature values are reproduced. Further details on the lattice setup are given in the Supplemental Material [14].
Naive measurements of G E (τ, T ) are highly susceptible to high-frequency fluctuations in the gauge fields and exhibit a poor signal-to-noise ratio. In quenched QCD, the multilevel algorithm [15] has been applied to overcome this problem. However, this algorithm is not applicable for QCD with dynamical fermions. To overcome the noise problem for our calculations with dynamical fermions we use the Symanzik-improved [16] gradient flow [17]. In quenched QCD it was demonstrated [18,19] that this approach is as effective as the multi-level algorithm for noise reduction, while also renormalizing G E nonperturbatively. By evolving the gauge fields in the fictitious flow time, τ F , as dictated by the force given by the gradient of the gauge action, the gradient flow smears the gauge fields over the radius √ 8τ F . Renormalization artifacts of the electric field operators E i due to finite lattice spacing a are highly suppressed for √ 8τ F > a. However, the flow radius should always be smaller than the relevant physical scales, implying the constraint √ 8τ F < τ < 1/(2T ). For G E it was found that the more strict criterion √ 8τ F /τ < 1/3 should be respected [18][19][20].
Results.-Since ρ E ∝ ω 3 for large ω, G E is a steeply falling function of τ . Therefore, it is convenient to normalize it by the leading-order perturbative result [7], with the Casimir factor, C F = 4/3, and the coupling constant, g, scaled out, that is, with G norm ≡ G LO E /(g 2 C F ). At small τ the lattice results will suffer from significant discretization effects. Furthermore, the distortions of the correlation functions due to gradient flow are the largest at small τ . The cutoff effects as well as the distortions due to gradient flow are also present in the free field theory. We can use the free theory result to estimate and partly correct for these effects. To reduce lattice cutoff effects as well as distortions due to gradient flow we perform tree-level improvement, meaning that we multiply the chromoelectric correlator by the ratio of the free correlator obtained in the continuum and  the one calculated on the lattice (in perturbation theory) with the given N τ at nonzero flow time: The details of calculating G norm τF (τ T, N τ ) can be found in Supplementary material [14].
The lattice chromo-electric correlators after tree-level improvement and normalized by G norm for the 96 3 × N τ lattices are shown in Fig. 1. We show the results for two different amounts of gradient flow, adjusted for each separation τ by fixing √ 8τ F /τ . In addition, we show the results from the coarsest lattices (64 3 ) as open symbols. We see that gradient flow is effective in reducing UV noise even for the largest lattice with N τ = 36. After treelevel improvement, the difference of G E /G norm obtained on the finest and the coarsest lattice is generally smaller than the statistical errors of the data obtained on the finest lattice. The flow time dependence is also quite small for τ T > 0.25 if the ratio of the flow radius √ 8τ F to τ is between 0.25 and 0.3. For τ T < 0.25 the amount of flow necessary to suppress discretization artifacts already comes close to the relevant physical scale of τ /3, leading to large distortions. For this reason the corresponding data points need to be omitted from the analysis.
Naively one expects that at high (but not extremely high) temperatures G E /G norm should not be different from unity since G LO E ≈ g 2 C F is a good approximation for G E and C F g 2 1. An interesting feature of the results shown in Fig. 1 is that the ratio G E /G norm has a much larger deviation from one than in the quenched case. In quenched QCD, G E /G norm reaches a value of about 4 at most [21]. This is due to the fact that the τ values in physical units (fm) accessible in full QCD are larger and, as we will see later, the value of κ in temperature units also turn out to be larger. As in quenched QCD, deviations from unity of G E /G norm are the largest at the lowest temperature, and become smaller as the temperature increases.
Next, we perform the continuum and flow-time-to-zero extrapolation of the chromoelectric correlator. First, we interpolate the correlators obtained on 64 3 × N τ lattices in τ for various values of τ F /τ 2 . From these interpolations we determine the correlator on coarser lattices for values of τ that are available for the finest 96 3 × N τ lattices and then perform continuum extrapolations for each τ T and τ F /τ 2 . As is apparent from Fig. 1, the cutoff effects are small except for small values of τ T . We perform the continuum extrapolation assuming that discretization errors go like (aT ) 2 ∼ 1/N 2 τ , which turns out to be capable of describing our data well, see Supplemental Material [14].
Finally, we perform the flow-time-to-zero extrapolation of the chromoelectric correlators. In the region a √ 8τ F τ we expect a linear τ F dependence as suggested by NLO perturbation theory [22]. And indeed, for 0.25 < √ 8τ F /τ < 0.3 a linear dependence seems to describe the data. Therefore, we use a linear extrapolation in τ F in this region to obtain the zero flow time limit, see Supplemental Material [14]. The continuum and zero flow time extrapolated results for the chromoelectric correlators are shown in Fig. 2. The extrapolations do not change the qualitative features of the correlation function but lead to a significant increase of the statistical errors.
With the continuum and flow-time-extrapolated data for the chromoelectric correlator we are in the position to estimate the heavy quark diffusion coefficient κ. To do so we need a parametrization of the spectral function that enters Eq. (2). Any parametrization of the spectral function should take into account its known behavior at small and large ω. For small ω the spectral function is solely determined by the heavy quark diffusion coefficient and has the form [7]: ρ E (ω, T ) ρ IR (ω, T ) = κω/(2T ) while at sufficiently large frequency the ω dependence of the spectral function should be described by perturbation theory due to asymptotic freedom in QCD. Moreover, thermal corrections to the spectral function are very small for ω T . Therefore, we assume that at large energies the spectral function is given by the LO or NLO perturbative T = 0 result up to a constant: The factor K accounts for the fact that the perturbative calculations may not be quantitatively reliable due to missing contributions from higher orders.
Perturbative calculations at NLO [7,23], classical simulations in effective three-dimensional theory [24] and strong coupling calculations [25] show that the spectral function is a smooth monotonically rising function of ω. Based on this, as well as the above considerations, we use the following two forms of the spectral function in our analysis that also have been used already in quenched QCD [21,26]: ρ max = max ρ IR (ω, T ), ρ UV (ω) and , which we refer to as the maximum (max) and the smooth maximum (smax) Ansätze, respectively. The latter is consistent with the perturbative NLO calculation [23] and OPE considerations [27] when it comes to the leading thermal correction at ω T . We also consider a third Ansatz for the spectral function, that is given by ρ IR (ω, T ) up to ω = ω IR , and by ρ UV for ω > ω UV , and for ω IR < ω < ω UV we interpolate with a power-law form ρ(ω, T ) = c ω p . The parameters c and p are chosen such that the spectral function is continuous at ω = ω IR and ω UV . This form of the spectral function has been used in Ref. [21]. Based on theoretical results we choose ω IR = T and ω UV = 2πT , see Supplemental Material [14].
Using the above three Ansätze for the spectral functions and the spectral representation of the chromoelectric correlator we fitted the continuum-and flow-timeextrapolated results treating κ and K as fit parameters and thus estimated the heavy quark diffusion coefficient. It turns out that the maximum Ansatz gives the largest value of κ, while the power-law form gives the smallest value. Using the LO or NLO form of ρ UV does not lead to significant change in the value of κ, meaning that the estimated values of κ are not too sensitive to the modeling of the high energy part of the spectral function.
Each model is fitted onto the same 1000 bootstrap samples of the double-extrapolated correlator data. We collect all results from all models in a single "distribution" for the fit parameter κ/T 3 . We determine a confidence interval by considering the median of this distribution, and then adding or subtracting the 34th percentiles on each side, which gives the lower and upper bounds of the interval. For better readability we quote the central value of the interval with the distance to the bounds as the uncertainty. We obtain: As already noted above, the shape of the correlation function for τ T > 0.25 does not seem to be significantly changed by the continuum and flow-time-tozero extrapolations. Therefore, we also performed the above analysis using the nonzero lattice spacing data at 1/a = 7.036 GeV and nonzero relative flow times √ 8τ F /τ = 0.3. We find that the estimated values of κ agree with the ones obtained from the continuumand flow-time-extrapolated data within errors. This is due to the fact that the systematic uncertainties associated with modeling of the spectral function are much larger than the effect of the continuum and zero-flowtime extrapolation. For this reason we also estimate the heavy quark diffusion coefficient at nonzero lattice spacing and flow time for the highest temperature resulting in κ( Conclusion.-We carried out first lattice QCD calculations of the heavy quark diffusion constant in 2+1 flavor QCD at leading order in the inverse heavy quark mass and in the phenomenologically relevant region of 195 MeV < T < 352 MeV. Our results for D s as function of T /T c are summarized in Fig. 3. Here we use T c = 180 MeV because the calculations are performed at m π 320 MeV, see Supplemental Material [14]. Our results are smaller than the quenched lattice QCD results [21,28]. At the lowest temperature our result agrees, within errors, with the strong coupling expecta-  tions from AdS/CFT [6,33]. At the highest temperature our result D s is compatible with the NLO perturbative prediction [5] within the uncertainties. In comparisons to some phenomenological determinations, namely, the Bayesian analysis [31] of heavy-ion collision data, and the ALICE Collaboration's model fits to their data [32], the lattice QCD results for D s are systematically smaller. On the other hand the the T -matrix approach on D s [29,30] seems to agree with the lattice results.
The present study can be extended in two different ways. Based on the previous lattice QCD studies of QCD equation of state [13], quark number susceptibilities [34,35], and static quark free energies [36] with light quark mass m l = m s /5 compared to those with nearly physical light quark mass m l m s /20, we expect the effect larger than physical quark mass on D s to be small in the temperature range considered in this study. However, such effects might be significant for temperatures closer to the QCD crossover temperature. Thus, we plan to extend the present calculations to smaller temperatures with physical values of the light quark masses. The heavy quark mass suppressed effects to the heavy quark diffusion coefficients are expected to be relatively small based on the calculations performed in quenched QCD [19,37]. We plan to estimate these corrections also in 2+1 flavor QCD.
All data from our calculations, presented in the figures of this paper, can be found in [38].
The computations in this work were performed using SIMULATeQCD [39,40]. Parts of the computations in this work were performed on the GPU cluster at Bielefeld University. We thank the Bielefeld HPC.NRW team for their support.
This material is based upon work supported by We thank the Institute for Nuclear Theory at the University of Washington for its kind hospitality and stimulating research environment during the completion of this work. INT is supported in part by the U.S. Department of Energy Grant No. DE-FG02-00ER41132.
We thank Guy D. Moore for valuable discussions.

Supplemental materials
In these supplemental materials we discuss important technical details on the calculation of the chromo-electric correlator and the heavy quark diffusion coefficient presented in the letter. In section I we discuss additional details about the lattice setup. In section II we discuss the calculation of the free chromo-electric correlator (leading-order result) on the lattice with gradient flow. In section III, we discuss the continuum and flow time extrapolations of the lattice results. Finally, in section IV, we discuss the reconstruction of the spectral function and the determination of the heavy quark momentum diffusion coefficient.

I. LATTICE SETUP
As mentioned in the paper, we performed lattice QCD calculations of the chromo-electric correlator using the HISQ action for quarks and the Symanzik-improved Lüscher-Weisz action for gluons. The calculations have been performed at the physical strange quark mass and light quark masses m l = m s /5 in a fixed scale approach on 96 3 ×N τ lattices for the bare gauge coupling β = 10/g 2 0 = 8.249, corresponding to a lattice spacing of 1/a = 7.036 GeV. The corresponding bare strange strange quark mass for this value of β in lattice units is am s = 0.01011 [12]. The temperature was varied by using N τ = 36, 32, 28, 25 and 20, corresponding to temperatures T = 195, 220, 251, 293 and 352 MeV, respectively. The lattice spacing and thus the temperature scale has been fixed using the r 1 -scale determined in Ref. [13] with the value r 1 = 0.3106 fm obtained in Ref. [41]. The value of the strange quark mass was obtained from the parametrization of the line of constant physics from Ref. [12]. The parameters of these calculations, including the statistics, are presented in Tab. I.
To control discretization effects and perform continuum extrapolations for each temperature (except the largest one), we performed calculations on two coarser lattices with lattice sizes 64 3 × N τ and N τ = 20, 22/24. Here, the temperature values were adjusted by varying the bare gauge coupling β. Again we used the r 1 -scale to set the temperature scale and lines of constant physics to obtain the quark masses from Ref. [13] and [12], respectively. In Tab. II we present the lattice details of the corresponding calculations including the temperatures, temporal lattice sizes N τ , quark masses and bare lattice gauge coupling, as well as the corresponding statistics.
The gauge configurations are generated using the RHMC algorithm. We save configurations every 10 trajectories and the acceptance rate is tuned to ≈ 80%. To estimate the auto-correlation time we consider the expectation value of the Polyakov loop evaluated on the flowed gauge configurations at the maximum flow time that enters the analysis, which is √ 8τ F T = 0.15. The integrated autocorrelation times turn out to be around 10-30 saved configurations. In each stream we bin the data into bins with length of one integrated autocorrelation time and use the approximately independent bins for the bootstrap analysis.
To facilitate the comparison of our lattice results on κ with available quenched QCD results and some model calculations, it is worth to show them as a function of temperature T in units of a (pseudo)critical transition temperature T c , chosen appropriately for each case. In the quenched case, T c is the critical temperature of the first-order deconfinement phase transition. In full QCD, this transition is only well defined at significantly larger light quark masses than what is used in this study. However, for the case of 2+1-flavor QCD with m l = m s /5 the light quark mass is sufficiently small to associate T c with the chiral crossover temperature. Typically the chiral crossover temperature is defined as the peak position in the disconnected chiral susceptibility. For the HISQ action with 2+1-flavor QCD and m l = m s /5 the disconnected chiral susceptibility has been calculated in Ref. [35] on 24 3 × 6 lattices. Using the corresponding data from Ref. [35] and performing fits with Gaussian and polynomial forms in the vicinity of the peak position, we obtain T c = 180(2) MeV, which we use to position the 2+1-flavor data in Fig. 3. The quoted error is statistical. The lower limit of the fit range was set to be T = 170 MeV, while the upper limit of the fit range was varied from T = 195 MeV to T = 210 MeV. The peak position does not depend on the fit range or on the form of the fit function (polynomial or Gaussian).

II. LEADING ORDER CORRELATION FUNCTIONS AT FINITE FLOW TIME IN LATTICE PERTURBATION THEORY
In the following, we will lay out the calculation of the leading order lattice perturbation theory correlation function G norm τF (τ T, N τ ) under flow that is used for tree level improvement [42]. Under gradient flow, the lattice gluon propagator takes the form where K (f ) is the flow kernel and D is the unflowed lattice propagator which is the inverse of the action kernel K (a) . The matrices depend on the discretization of gradient flow and action respectively. For Zeuthen flow and the Lüscher- Weisz gauge action, the kernels can be found in Ref. [16]. Using the flowed propagator, leading order perturbation theory expectation values of some gluonic operator M under flow can be written as where M AB µν is an observable-dependent matrix, and . . . 0 means expectation value in the free field theory. The symbol p is to be understood as a shorthand notation for T p0 d 3 p/(2π) 3 . By taking the lattice expression of the color-electric field strength G E and expanding the link variables U µ (x) = exp(iag 0 A µ (x)) in the Lie algebra fields, we find the operator matrix Here, with the lattice momentump µ = (2/a) sin(p µ a/2). Using this, Eq. (5) is integrated numerically for each τ and τ F to obtain the leading-order result on the lattice, G LO E,τF (τ T, N τ ). The accuracy of the integration is tested by making use of the fact that for Wilson action and Wilson flow, an analytical solution to the integral exists [18]. Numerical integration for the Wilson-Wilson case agrees with the analytical result up to 10 −8 which is much smaller than the statistical uncertainty, so the integration procedure is well suited. The normalized lattice correlator is then defined as G norm

III. CONTINUUM AND FLOW TIME EXTRAPOLATIONS
As in the main text in the following discussion we will normalize the chromo-electric correlator by the free theory (LO) result with the Casimir factor C F = 4/3 and the coupling constant g 2 scaled out [7]: .  To perform the continuum extrapolation we first use spline interpolations of the 64 3 × N τ (N τ = 20, 22/24) lattice data for G E in τ T at each temperature. From the spline interpolations we determine G E at the values of τ T available on the finest (96 3 × N τ ) lattices for the same temperature. In this way we can study the lattice spacing dependence of G E /G norm at each temperature in terms of aT = 1/N τ . This is shown in Fig. 4 for √ 8τ F /τ = 0.3 as an example. In general, we see that the lattice spacing (N τ ) dependence of G E /G norm is mild.
Unfortunately, even with gradient flow the statistical errors on the data are too large to perform naive individual continuum extrapolations for each distance τ T , as these will suffer from overfitting of the noise. However, from our perturbative lattice calculations and from previous quenched results we expect the continuum to be approached from below, and we expect the slope of the extrapolation to vanish for very large τ T , with a leading decay proportional to 1/(τ T ) 2 . Therefore, we constrain the slope of the extrapolation by performing a combined extrapolation of all τ T at fixed . This is shown in Fig. 4 for two temperatures at √ 8τ F /τ = 0.3. For all temperatures and all relative flow times this fit ansatz yields a χ 2 /d.o.f ∈ [0.8, 1.9], meaning that the data is described well and that this constraint is not too restrictive given the statistical precision. To estimate the error on the continuum extrapolated results we use the bootstrap method.
Next we examine the flow time dependence of the continuum extrapolated result for the chromo-electric correlator. The flow time dependence of the chromo-electric correlator is shown in Fig. 5. As mentioned in the main text, in the window a √ 8τ F τ /3 we expect a linear τ F -dependence suggested by the NLO perturbative calculation [22]. In the continuum extrapolated data we find that the upper bound for the flow time as suggested by perturbation theory is a bit too generous (meaning that nonlinear behavior starts already at slightly smaller flow times), which is why we restrict the upper limit to √ 8τ F /τ < 0.3. The lower bound is also chosen more conservatively, in order to be sure that discretization artifacts are well-suppressed and in order to obtain a good enough signal-to-ratio. For these reasons we perform the flow extrapolations in the range 0.25 < √ 8τ F /τ < 0.3 for τ T ≥ 0.25. From our perturbative results, quenched studies and the nonzero lattice spacing data of this study, we expect the flow time limit to be approached from below, and we also expect the slope of the extrapolation to be roughly independent of τ T . For these reasons, we again perform a combined extrapolation of all τ T with the Ansatz f (τ F ; τ T ) = τ T G τF=0 τ T + m τ F . Given our statistical precision in the continuum, this Ansatz can describe the data well for all temperatures.
We note that the results on G E evaluated for different flow times are strongly correlated. Therefore, for the zero flow time extrapolation we restrict ourselves to use only three τ F values as indicated in Fig. 5. Adding data on G E at more τ F values will not change the result of extrapolation. We note that here we use uncorrelated fits because given our statistics it is impossible to reliably estimate the correlation matrix. The errors on the extrapolated values are again estimated using the bootstrap procedure.
We note that we also performed a second analysis for the continuum and flow-time-to-zero extrapolation, where the continuum and flow time extrapolations are performed independently at each τ T value. In this analysis, to estimate the continuum limit we perform a constant fit of the G E /G norm data for τ T > 0.33, while for smaller τ T values we perform a fit with a linear Ansatz in 1/N 2 τ . As explained above, we expect the flow time limit to be approached from below, which is why we constrain the slope of the extrapolation for each τ T to be smaller or equal zero. The results for the continuum and zero flow time extrapolated correlator and κ/T 3 from this analysis agree with the above approach well within errors.

IV. SPECTRAL RECONSTRUCTION UV part of the spectral function models
In order to obtain the heavy quark diffusion coefficient κ, an Ansatz for the spectral function ρ E is needed. The behavior of the spectral function at large energies, denoted by ρ UV , is an essential input for all Ansätze of the spectral function. We use the perturbative LO and NLO results at T = 0 for ρ UV , multiplied by a free parameter K, which accounts for higher-order corrections (K-factor). To evaluate ρ UV , we need to fix the renormalization scale at which the coupling constant g 2 is defined. At LO, we use the naive choice for the scale µ = max(µ T , ω), where µ T is the thermal scale inferred from the high temperature three dimensional effective theory scale [43] with N c being the number of colors and N f being the number of light quark flavors. In the context of the chromoelectric spectral function, the use of this scale was suggested in Ref. [23]. For the case of interest, N c = 3 and N f = 3, µ T 9.08222T . The logarithmic term in the NLO correction sets the natural scale when the NLO result is used. In this case we use the "optimal" scale that sets the NLO correction to zero for large ω [23]: For N c = 3 and N f = 3, µ opt 14.7427ω is large for all relevant values of ω, so there is no need to assume that the renormalization scale needs to be temperature dependent. We use the 5-loop running coupling constant in our analysis as implemented in the "AlphasLam" routine of the RunDec package [44,45] and the value of Λ MeV [46][47][48][49][50][51][52][53][54].

Fit results
All of spectral function models considered in this study exhibit two fit parameters: the momentum diffusion coefficient κ/T 3 , and the K-factor that multiplies the UV part. The spectral functions corresponding to the fit results for each model are shown in Fig. 6 for the lowest and the highest temperatures (the intermediate temperatures look similar). The fits work well as can been seen in Fig. 7, where we show the fitted model correlators normalized to the data for different models.
The values of the K-factor are shown in Fig. 8 as a function of T /T c for the smax Ansatz as an example; the results for the other models are again similar. The value of the K-factor is close to one when the LO result is used for ρ UV , while it is around two when using the NLO result. The situation is somewhat similar to the one in quenched QCD. The analysis of Refs. [18,26] used the LO result for ρ UV and found K ∼ 1, while the analysis of Refs. [19,21] that used the NLO result for the UV part of the spectral finds K 1.7 at similar values of T /T c .
We also see from Fig. 6 that there is a difference between the LO and NLO forms at large ω, but at the same time this does not really seem to effect the slope of ρ E /ω in the IR region, meaning that the determination of κ/T 3 is not affected much by the modeling of the UV part of spectral function, in particular by the value of K. The different values of K for the LO and NLO forms of the spectral functions are largely due to the difference in the choice of the renormalization scale µ. For larger values of µ the K factor for LO form of ρ UV would be very different from one.
The κ/T 3 values obtained from all different fits are shown in Fig. 9 at the lowest and highest temperature. The temperature dependence of κ/T 3 is shown in Fig. 10 and Tab. III. We see that κ/T 3 decreases with increasing temperature as also suggest by the weak coupling calculations. We also see that the central values of κ/T 3 obtained in 2+1-flavor QCD are significantly larger compared to the ones in quenched QCD, see Refs. [18,19,21,26,28,55]. We note that quenched results on κ/T 3 from various lattice groups agree well (cf., Fig. 3 in Ref. [28]).

Details for power-law model
In the main text we introduced the power-law Ansatz, where ρ E (ω, T ) is equal to ρ IR (ω, T ) for ω < ω IR , and is equal to ρ UV (ω) for ω > ω UV . We need some prior knowledge to set ω IR and ω UV . From the NLO calculations we know that the thermal corrections to the spectral function are very small for ω > 2πT [23]. Classical calculations in the effective three dimensional theory suggest that the spectral function ρ E is approximately linear in ω for ω < 0.8T [24]. These calculations are expected to describe the non-perturbative physics associated with soft chromo-magnetic fields at high temperatures. The strong coupling calculations based on AdS/CFT on the other hand show that the spectral function is approximately linear up to ω πT [25]. Therefore, for the temperature range of interest, ω IR should not be much smaller than one. We choose ω IR = T in our analysis. The corresponding fits also work well as one can see from Fig. 7. The K factors for these fits are very similar to the ones obtained for max and smax forms.
The spectral functions corresponding to the power law Ansatz are also shown in Fig. 6 at the highest and lowest temperatures. As one can see from the figure, this Ansatz leads to smaller y-intercepts of ρ E /ω, and thus smaller values of κ/T 3 . Choosing larger values of ω IR would result in κ/T 3 that is closer to the ones obtain from the max and smax Ansätze. As one can see from Fig. 6, for the max and smax Ansätze, the linear behavior of the spectral function extends to ω ≥ πT . Therefore, the max and smax Ansätze correspond to spectral functions that are more similar to the spectral function in the strongly coupled theory, while the power law Ansatz with ω IR = T is closer to the spectral functions obtained in the weak coupling approach. Thus the max Ansatz and the power-law Ansatz for the spectral function could be interpreted as incorporating features from the strong and weak coupling calculations, respectively.

Systematics at nonzero lattice spacing and flow time
In the main text we saw that the overall shape of the chromo-electric correlator does not change much after the continuum and zero flow time extrapolation (cf., Fig. 1 and 2). The effects of nonzero lattice spacing and nonzero flow time mostly affect the chromo-electric correlation function at small τ T , which are therefore not used in this analysis. Since we cannot perform a continuum extrapolation for the highest temperature (T = 352 MeV), due to a lack of finer lattices, the question arises how much of a systematic error is introduced when performing fits on non-extrapolated correlator data. This is further incentivised by considering a generalization of the spectral representation in Eq. 2 to nonzero lattice spacing, which implies that the corresponding spectral function has support only up to some maximal energy ω max ∼ 1/a. For local (point-like) meson operators this has been shown in Ref. [56]. If one considers extended instead of point-like meson operators, the ω max of the corresponding spectral function will decrease [57]. However, the spectral functions of the extended meson at ω significantly smaller than the cutoff will be the same as for the correlator of a local meson [57], if the extent of the meson operators is of the order of few lattice spacings.
Correlators calculated with the gradient flow can be considered as correlators of extended operators as the gradient flow smears the field within the radius of √ 8τ F . Based on these considerations we expect that the non-zero lattice spacing and flow time will mostly affect the large ω behavior of ρ E (ω, T ) and only have a small effect on its low ω behavior and thus the value of κ/T 3 . For quenched QCD this has been shown in Ref. [19]. Therefore, we performed the above spectral analysis also for the chromo-electric correlator at a −1 = 7.036 GeV and flow time √ 8τ F /τ = 0.3. As expected, the fitted spectral function is different in the UV region, which manifest itself in the smaller value of the K factor as one can see from Fig. 8. On the other hand the values of κ/T 3 obtained from the analysis of G E at a −1 = 7.036 GeV and √ 8τ F /τ = 0.3 should not change much. As shown in Fig. 10 and Tab. III, the κ/T 3 values obtained from the continuum and zero flow time extrapolated chromo-electric correlator agree well with the ones obtained at a −1 = 7.036 GeV and √ 8τ F /τ = 0.3 within errors, meaning that the statistical and the systematic uncertainties in κ/T 3 are much larger than the effects of non-zero lattice spacing and gradient flow. For this reason we extend our lattice estimate of κ/T 3 also to T = 352 MeV, which is shown in Fig. 10 and Tab. III. ρ model = ρ "smax"    Fig. 10. The second column shows the result obtained from the fits of the continuum and zero flow time extrapolated results of GE. The third column shows the results for κ/T 3 obtained from the fits of the chromoelectric correlator calculated with a −1 = 7.036 GeV at relative flow radius √ 8τF/τ = 0.3.