Lattice QCD constraints on the heavy quark diffusion coefficient

We report progress towards computing the heavy quark momentum diffusion coefficient from the correlator of two chromo-electric fields attached to a Polyakov loop in pure SU(3) gauge theory. Using a multilevel algorithm and tree-level improvement, we study the behavior of the diffusion coefficient as a function of temperature in the wide range $1.1<T / T_c<10^4$ in order to compare it to perturbative expansions at high temperature. We find that within errors the lattice results are remarkably compatible with the next-to-leading order perturbative result.


I. INTRODUCTION
The matter produced in heavy ion collisions can be described as a nearly ideal fluid, see Ref. [1] for a recent review. Because of the high energy density, the created matter is deconfined and can be characterized as a strongly coupled quark-gluon plasma (sQGP) [2,3]. One recently realized interesting feature of the quark gluon plasma is the fact that heavy quarks participate in the collective behavior, see Ref. [4] for a recent review. This is interesting for the following reason. The relaxation time of heavy quarks is expected to be ∼ (M/T )t light rel , with M being the heavy quark mass, T being the temperature, and t light rel being the relaxation time of the bulk (light) degrees of freedom in sQGP. The lifetime of the hot medium created in heavy ion collisions is about 5−10 fm. Since the collectivity in the heavy quark sector implies that the relaxation time of the heavy quark is much shorter than the lifetime of the medium despite the enhancement factor of M/T , this in turn means that the relaxation time of the bulk degrees of freedom is very short, thus further corroborating the strongly coupled nature of the matter produced in heavy ion collisions.
Because the relaxation time of heavy quarks is much larger than the relaxation time of light degrees of freedom the dynamics of heavy quarks can be understood in terms of Langevin equations [5]. The drag coefficient η and heavy quark momentum diffusion coefficient κ that enter into the Langevin equations describe the interaction of the heavy quarks with the medium and are connected by the Einstein relation η = κ/(2M T ) in thermal equilibrium. The heavy quark diffusion coefficient has been calculated in perturbation theory at leading order (LO) [5,6] as well as at next-to-leading order (NLO) [7]. * nora.brambilla@ph.tum.de † viljami.leino@tum.de ‡ petreczk@bnl.gov § antonio.vairo@ph.tum.de The NLO correction is very large, thus questioning the validity of the perturbative expansion. Analytic calculations for strong coupling are available only for supersymmetric Yang-Mills theories [8,9]. Therefore, lattice QCD calculations for the heavy quark diffusion coefficient are needed.
It is well known, however, that lattice calculations of the transport coefficients are very difficult. To obtain the transport coefficients one has to reconstruct the spectral functions from the appropriate Euclidean time correlation functions. At low energies, ω, the spectral function has a peak, called the transport peak, and the width of the transport peak defines the transport coefficient. Thus, one needs a reliable determination of the width of the transport peak in order to obtain the transport coefficient from lattice QCD calculations, which is difficult [10,11]. In the case of heavy quarks, this is even more challenging because the width of the transport peak is inversely proportional to the heavy quark mass. Moreover, Euclidean time correlators are rather insensitive to small widths [11][12][13][14][15]. Recently the problem of heavy quark diffusion has also been studied out of equilibrium with real-time lattice simulations in Refs. [16,17]. Moreover, the heavy quark momentum diffusion coefficient is a crucial parameter entering the evolution equations describing the out-of-equilibrium dynamics of heavy quarkonium in sQGP [18][19][20].
The above difficulty in the determination of the heavy quark diffusion coefficient can be circumvented by using an effective field theory approach. Namely, by integrating out the heavy quark fields one can relate the heavy quark diffusion coefficient to the correlator of the chromoelectric field strength [21]. The corresponding spectral function does not have a transport peak and the small ω behavior is smoothly connected to the UV behavior of the spectral function [21]. The heavy quark diffusion coefficient is given by the intercept of the spectral function at ω = 0 and no determination of the width of the transport peak is needed. Lattice calculations of κ along these lines have been carried out in the SU(3) gauge theory in the deconfined phase, i.e. for purely gluonic plasma [22][23][24][25][26]. The correlator of the chromo-electric field strength is very noisy making the lattice calculations extremely challenging. To deal with this problem it is mandatory to use noise reducing techniques such as the multi-level algorithm by Lüscher and Weisz [27]. This algorithm is based on the locality of the action and therefore is only available for the pure gauge theory. This is the reason why the calculations of the heavy quark diffusion coefficient are performed in the SU(3) gauge theory. Another challenge in the determination of the heavy quark diffusion coefficient is the reconstruction of the spectral function from the Euclidean time correlation function. The above lattice studies used a simple parameterization of the spectral function to extract κ. One has to explore the sensitivity of the results on the parameterization of the spectral function. More generally, one has to understand to what extent the Euclidean time correlation function of the chromo-electric field strength is sensitive to the small ω behavior of the corresponding spectral function.
At sufficiently high temperatures the perturbative calculations of the heavy quark diffusion coefficient should be adequate. This suggests that κ/T 3 should decrease from large values at temperatures close to the transition temperature to smaller values when the temperature is increasing. It would be interesting to see if contacts between the lattice and the perturbative calculations can be made for the heavy quark diffusion coefficient as it has already be done for the equation of state [28], quark number susceptibilities [29,30] and static correlation functions [31][32][33]. If such contacts can be established these would validate the methodology used in the lattice extraction of κ. Previous lattice studies focused on a narrow temperature region [24] or only considered a single value of the temperature [25]. In Ref. [24] no significant temperature dependence of κ/T 3 was found. Large temperatures are needed in the lattice studies to establish the temperature dependence of κ/T 3 . The temperature dependence of κ/T 3 is also important for phenomenology as with a constant value of κ/T 3 it is impossible to explain simultaneously the elliptic flow parameter, v 2 , for heavy quarks and the nuclear modification factor [4]. Furthermore, the spectral function of the chromo-electric field strength correlator is known at NLO [34]. Using this NLO result at high ω one can constrain the functional form of the spectral function used in the analysis of the lattice correlator.
The aim of this paper is to study the correlator of the chromo-electric field strength in a wide temperature range in order to make contact with weak coupling calculations of the Euclidean correlation function up to NLO in the spectral function, and also to constrain the temperature dependence of κ.
The rest of the paper is organized as follows. In the next section we go trough the procedure of calculating the Euclidean correlator of the chromo-electric field strength on the lattice. The spectral function of the chromoelectric correlator and its relation to κ is discussed in section III. There we also review the perturbative results for this spectral function. The short time behavior of the chromo-electric correlator and its proper normalization is clarified in section IV. In section V, we discuss how to model the spectral functions of the chromo-electric correlator and to extract the value of κ from the lattice results. Finally, section VI contains our conclusions.

II. LATTICE RESULTS FOR THE CHROMO-ELECTRIC CORRELATOR
For a heavy quark of mass M πT , the heavy quark effective theory (HQEFT) provides a method of calculating the heavy quark diffusion coefficient in the heavy quark limit by relating it to a chromo-electric correlator in Euclidean time [9,21]: where T is the temperature, U (τ 1 , τ 2 ) is the temporal Wilson line between τ 1 and τ 2 , and the chromo-electric field, in which the coupling has been absorbed E i ≡ gE i , is discretized on the lattice as [21]: This discretization is expected to be least sensitive to ultraviolet effects [21].
To calculate the discretized chromo-electric correlator defined above on the lattice we use the standard Wilson gauge action and the multilevel algorithm [27]. We consider N 3 s × N t lattices and vary the temperature in a wide range T = 1.1 T c − 10 4 T c by varying the lattice gauge coupling β = 6/g 2 0 . Here T c is the deconfinement phase transition temperature. We use N t = 12, 16, 20 and 24 at each temperature to check for lattice spacing effects and perform the continuum extrapolation. In this study we use N s = 48, except for N t = 12 lattices, where multiple spatial volumes are used to check for finite volume effects.
To set the temperature scale as well as the lattice spacing we use the gradient flow parameter t 0 [35] and the value T c √ t 0 = 0.2489(14) [36]. We use the result of Ref. [36] to relate the temperature scale or the lattice spacing to β. The parameters of the lattice calculations including the statistics are given in Table I. In the simulations with the multilevel algorithm we divide the lattice into four sub-lattices and update each sub-lattice 2000 times to evaluate the chromo-electric correlator on a single gauge configuration. We use the simulation program developed in a prior study [24].
In order to obtain the heavy quark diffusion coefficient, the lattice chromo-electric correlator needs to be renormalized and then extrapolated to the continuum. 1 The renormalization coefficient Z E (β) ≡ Z E (g 2 0 ) of the chromo-electric correlator in case of the Wilson gauge action has been calculated at 1-loop [37]: We will use this 1-loop correction in the present study. However, we expect that the 1-loop result for Z E is not precise enough. As it will be clear from the results of the lattice calculations this is indeed the case. The perturbative error in Z E (β) affects both its absolute value for fixed β and its β-dependence. For the continuum extrapolation it is important to estimate the uncertainty in the β-dependence of the renormalization constant. The error in the absolute value of Z E could be corrected after the continuum extrapolation is done by introducing an additional multiplicative factor. We will postpone the discussion of this multiplicative factor to section IV. To estimate the error in the β dependence of Z E we consider the tadpole improved result for Z E , namely Z tad E = 1/u 0 , with u 0 being the plaquette expectation value [24]. The difference in the β dependence of Z tad E and Z 1−loop E can be used as an estimate of the error of the β-dependence of Z E . Therefore, at each temperature we consider the variation in Z 1−loop E · u 0 in the β range that corresponds to N t = 12 − 24 as an estimate of the systematic errors in Z E for bare gauge couplings in that range.
The chromo-electric correlator decays rapidly with increasing τ . This feature can be understood from the leading order (tree-level) result [21]: where C F = 4/3 is the Casimir of the fundamental representation of SU(3). In Fig. 1 we show Z E G E /G norm E for different temperatures calculated on the largest, 48 3 × 24 lattice. We see a significant temperature dependence in this ratio. Also shown in the figure are the numerical results for the lowest temperature, T = 1.1 T c calculated for different N t . As one can see from the figure, the cutoff (N t ) dependence is significant even for relatively large values of τ T . We expect that the cutoff dependence increases with decreasing τ T , except when τ is of the order of the lattice spacing because the cutoff dependence of We see that our lattice data follow this expectation for τ T > 0.2. This observation is important for estimating the reliability of the continuum extrapolations. A similar N t dependence is observed at other temperatures. In order to reduce discretization errors we turn to a tree-level improvement procedure [38,39], where the leading order results in the continuum (4) and the lattice perturbation theory are matched. The LO lattice perturbation theory gives [23]: for all temperatures when tree-level improvement is used.
As an example, we show this reduction at the bottom of Fig. 2 for the lowest temperature T = 1.1 T c . A similar reduction in the N t dependence is seen at other temperatures. Due to its impact, we will use the tree level improvement for the rest of this paper and, therefore, unless otherwise indicated, drop the overline from τ and the superscript imp from G imp E . The normalized chromo-electric correlator shown in Fig. 2 has a significant τ -dependence. We conclude that the LO perturbative result does not capture the key features of the chromo-electric correlator. Only at the highest temperature, T = 10 4 T c , the τ -dependence of the correlator is well described by the leading order result. One may wonder whether the observed behaviour of the normalized chromo-electric correlator is due to thermal effects that are not present at leading order, like the physics of the heavy quark transport or are due to higher order effects at zero temperature. In order to answer this ques- The chromo-electric field correlator from Eq. (1) normalized with Eq. (4) and tree-level improved with (5). Top: All measured temperatures for the biggest lattice size. Bottom: All lattice sizes for the smallest temperature.  tion we show our lattice results in Fig. 2 as function of τ in physical units rather than of τ T . This figure shows that the ratio of the chromo-electric correlator to the free theory result is largely temperature independent implying that the chromo-electric correlator is dominated by the vacuum part of the spectral function. Even more important it becomes to quantify the temperature dependence of the chromo-electric correlator. This can be done by considering the following ratio of the normalized correlator at fixed value of β but at two temperatures corresponding to temporal extent N t and 2N t . Lattice artifacts are canceled out in the double ratio: We present our results for this quantity in Fig. 4 and observe that in the small τ T region there is no temperature dependence. Instead the thermal effects become more noticeable as τ T is increased, being highest at the largest τ T . Also, these thermal effects become more prominent as the temperature is lowered towards the transition temperature. In particular, we see a near negligible temperature dependence in the T = 10 4 T c data. As we will see in the next sections, these features of R 2 (N t ) can be naturally explained by a temperature dependent κ. In any case thermal effects in the chromo-electric correlator, which also encode the value of κ, are small, at the level of few percent. This fact implies that extracting κ from lattice determinations of the chromo-electric correlator is challenging. Before extracting the heavy quark diffusion coefficient we need to address finite volume effects and perform the continuum extrapolation of Z E G E . Most of our calculations have been performed using N s = 48. To check for finite volume effects for N t = 12 we have performed calculations using spatial sizes N 3 s = 24 3 , 32 3 , 48 3 at two temperatures, T = 1.5 T c and T = 10 T c . The smallest spatial volume here corresponds to the aspect ratio N s /N t = 2. The detailed study of finite volume effects is discussed in Appendix A. We find that the finite volume effects are small, considerably smaller than other sources of errors down to the aspect ratio N s /N t = 2. Therefore, at the current level of precision using a N s = 48 lattice is sufficient even for N t = 24.
Next, we perform the continuum extrapolations of Z E G E . The systematic errors in the renormalization constant estimated above are combined with the statistical errors of the chromo-electric correlator before performing the continuum extrapolation. In the interval 0.1 ≤ τ T ≤ 0.45 we have a sufficient number of data points to perform the continuum extrapolations. We first interpolate the data for each N t in τ T using ninth order polynomials to estimate Z E G E at common τ T values. We perform linear extrapolations in 1/N 2 t = (aT ) 2 of Z E G E at these τ T values using lattices with N t = 16, 20 and 24. As an example, we show the continuum extrapolation for selected values of τ T in Fig. 5. One can see that the N t = 12 data does not lie in the 1/N 2 t scaling region. Therefore, we also perform extrapolations to N t = 12 data with a (aT ) 4 term included. The difference between these continuum extrapolations is used as an estimate of the systematic error of the continuum result. The slope of the a 2 dependence is increasing with decreasing τ T as can be seen from Fig. 5. This is expected, the cutoff effects are larger at smaller τ T . However, at the smallest value τ T = 0.1 the slope of the a 2 dependence becomes smaller again contrary to the expectations. We take this as an indication that the cutoff effects in this region cannot be described by a simple a 2 or a 2 + a 4 . As shown in Appendix A the slope of the a 2 dependence increases monotonically only till τ T ≥ 0.175. Therefore, we consider the continuum extrapolation to be reliable only for τ T ≥ 0.175. For an additional cross-check, we also perform the continuum extrapolation of the lattice data without tree-level improvement. This is discussed in Appendix A, where further details of the continuum extrapolations can be found.
The continuum extrapolated chromo-electric correlator normalized by G norm E is shown in Fig. 6 for all temperatures as function of τ T . The continuum extrapolated results share the general features of the tree-level improved results at non-zero lattice spacing in terms of τ and temperature dependence. In particular, we see a strong dependence on τ T , except for the highest temperature, indicating that the leading order result does not capture the τ T dependence of G E . We will try to understand these features of the correlator in the next sections.

III. SPECTRAL FUNCTIONS AND DIFFUSION COEFFICIENT IN PERTURBATION THEORY
In order to determine the heavy quark diffusion coefficient κ from the chromo-electric correlator, G E one has to use the relation between this correlator and the spectral where The heavy quark diffusion coefficient is determined in terms of ρ through the Kubo formula [40] κ ≡ lim ω→0 2T ρ(ω, T ) ω .
At leading order of perturbation theory the spectral function is given by [21]: where the coupling has been evaluated at the scale µ ω . We use the 5-loop running coupling constant in this work [41]. At LO the scale µ ω is arbitrary. A natural choice is µ simple ω = max(ω, πT ) [25]. The LO spectral function (11) gives κ = 0.
At NLO the perturbative calculation of ρ(ω, T ) needs Hard-Thermal-Loop (HTL) resummation for ω < ∼ m E , with m E being the LO Debye mass: m E = N c /3gT in the pure gauge theory. The full NLO result of ρ(ω, T ) has been calculated in [34]. The NLO spectral function provides the LO non-vanishing result for κ: where For ω > ∼ T there is no need for resummation when calculating the spectral function at NLO; the naive (non-resummed) NLO result for ρ(ω, T ) in the pure gauge case reads where n B (q) = (exp(q/T )−1) −1 is the Bose-Einstein distribution, P takes the principal value, and g 2 ≡ g 2 (µ ω ).
The first line of (13) gives the NLO T = 0 contribution and the subsequent lines carry the thermal effects. For the NLO ρ(ω, T ), µ ω may be set such that the NLO T = 0 contribution vanishes [34]: and the T = 0 part of (13) reduces to (11). This is a convenient choice of scale for ω T . For ω ∼ T or smaller a convenient choice of scale was proposed in Ref. [42] ln(µ ω ) = ln(4πT in the pure gauge case. We switch between these two scales when they become equal at ω 0.8903 T [34]. The heavy quark diffusion coefficient has been calculated at NLO and the result reads [7]: The NLO result for κ cannot be replicated from currently known spectral functions as that would require ρ(ω, T ) to be available at NNLO, which it is not. Both the LO and NLO result for κ are obtained under the weak coupling assumption m E T . This condition, however, is not satisfied for most of the temperatures of interest. As a consequence, one obtains an unphysical behavior at LO, i.e., that κ becomes negative for T < 10 3 T c .
One can also calculate κ using the kinetic theory. The corresponding expression reads [7,43]: If we do not expand in the temporal gluon self-energy, Π 00 (p), which is formally of order g 2 , the above expression contains higher order contributions to κ as well. Therefore, the above expression can be considered as the resummed leading order result. The temporal gluon selfenergy depends on the gauge choice. For small momenta it can be expanded as The first two terms in this expansion are gauge independent. We can take either the first term or the first and second terms in the above expression and evaluate the integral in Eq. (17) numerically. Only keeping the first term in the above expression for Π 00 already leads to a positive result, while keeping also the second term leads to an enhancement of the κ value. We present all the different perturbative results for κ as a function of temperature in Fig. 7. The scale of the coupling is the one defined in Eq. (15). At the highest temperature, T = 10 4 T c , considered in this study we expect that the NLO result can provide some guidance on the properties of the spectral function and on the τ dependence of the chromo-electric correlator. Therefore, in Fig. 8 we show different versions of the NLO spectral function, including the zero temperature one. The full NLO spectral function can be described well by the simple κ LO ω/(2T ) form for ω < 0.02 T , while it approximately agrees with the T = 0 result for ω > 2 T . The full NLO result and the naive (unresummed) NLO result agree for ω > 0.6 T . At small ω the naive NLO result is logarithmically divergent. This divergence cancels against contributions coming from the scale m E in the resummed expression. We can model the spectral function  by smoothly matching the κ LO ω/(2T ) behaviour at small ω with the zero temperature spectral function at large ω. We call this the perturbative step form. It is also shown in Fig. 8 by the blue dotted line. Using the NLO spectral function we can calculate the corresponding chromoelectric correlator, which is shown in Fig. 9. The width of the band corresponds to the variation of the scale by a factor two around the value given by Eq. (14). The scale variation appears to be very small. We also calculate the chromo-electric correlator corresponding to the perturbative step form. The resulting correlator is indistinguishable from the one obtained using the NLO spectral function. This means that the additional structures in the spectral function in the region 0.02 < ω/T < 0.6 play no significant role when it comes to the correlator. We have also considered a perturbative step model using κ NLO . While using the NLO result for κ significantly enhances the spectral function in the low ω region it only leads to a 0.2% enhancement of the chromo-electric correlator compared to the one obtained using κ LO . Thus, the correlator is not sensitive to the small ω part of the spectral function at the highest temperature. At lower temperatures κ gets larger and the contribution of the low ω part of the spectral functions is more prominent. Therefore, it is at lower temperatures that the value of κ can be constrained by accurate calculations of the chromo-electric correlator.
While at T = 10 4 T c , one may expect the resummed NLO result to provide an adequate description of the spectral function, this is not expected at lower temperatures, because, as pointed out above, numerically m E > T . In particular, for T < 10 3 T c the resummed spectral function turns negative at some point in the region ω < T , thus implying that the resummed perturbative result is not applicable in this ω range. In section V we will discuss the implications of this finding.
In Fig. 9 we also show the continuum limit of the chromo-electric correlator at the high temperature T = 10 4 T c for comparison. The continuum extrapolated lattice result of the chromo-electric correlator has the same shape as the NLO calculation. We note, however, that our continuum data differs from the perturbative curve by a factor 1.2, which indicates that the renormalization constant is not accurate. (which is the same up to a multiplicative factor) does not capture the τ dependence of the chromo-electric correlator. Therefore, as an alternative normalization we consider a correlator obtained from Eq. (9) using the zero temperature NLO result for the spectral function with running coupling constant evaluated at scale µ ω given by Eqs. (14) and (15). We label the corresponding correlator as G NLO+ are shown in Fig. 10. We see that this ratio increases less with increasing τ and there is also some indication of appearance of a plateau at small τ T . This indicates that G NLO+ E captures the τ dependence of the chromoelectric correlator obtained on the lattice much better. However, even at the smallest τ the ratio G E /G NLO+ E is different from one. This is most likely due to the fact that the 1-loop result is not accurate for Z E . As shown in the previous section, even for the highest temperature, T = 10 4 T c , the NLO result is lower by a factor 1.2, as seen in Fig. 9, although its τ -dependence agrees well with the continuum extrapolated lattice data. Therefore, we introduce an additional normalization factor, C N by normalizing the ration G E /G NLO+ E to one at τ T = 0.19. To check the uncertainty of C N due to the choice of the normalization point, we also consider τ T = 0.175 as possible normalization point. Furthermore, we vary the scale µ ω by a factor of two around the optimal value when evaluating C N . The numerical values of C N are shown in Tab. II for different temperatures. The dependence on the normalization point is shown in the systematic error and is of the same order as the scale dependence. The additional normalization constant C N decreases with increasing temperature. This is due to the fact that the β range used in the evaluation of the lattice correlator is increasing with increasing temperature and the 1-loop result is more reliable at large β values. We will normalize G E /G NLO+ E with C N given in Tab. II before comparing with the model spectral functions used for the extraction of κ.

V. MODELING THE SPECTRAL FUNCTION AND DETERMINATION OF κ
To obtain the heavy quark diffusion coefficient from the continuum extrapolated lattice results we need to assume some model for the spectral function. We will use the NLO results on the spectral function as well as κ to guide us in this process. We also need to consider how sensitive the Euclidean time chromo-electric correlator is to the spectral function in different ω regions. From the previous sections it is clear that G E is dominated by the large ω part of the spectral function and thermal effects in the spectral function contribute at the level of few percent to the correlator.
It is reasonable to assume that at large enough ω perturbation theory is reliable even if the condition m E T is not satisfied. This is because for large ω HTL resummation is not important, as will be detailed later. Certainly at zero temperature the perturbative calculation of ρ(ω, T ) is reliable for ω Λ QCD . Therefore, we assume that for ω > ω UV the spectral function is given by ρ UV (ω, T ), which is calculated perturbatively. On the other hand for sufficiently small ω the spectral function is given by and we can assume that ρ(ω, T ) = ρ IR (ω, T ) for ω < ω IR . In the region ω IR < ω < ω UV the form of the spectral function is not known, in general, and this lack of knowledge will generate an uncertainty in the determination of κ. We consider two possible forms of the spectral functions that are continuous and are based on simple interpolations between the small ω (IR) region and large ω (UV) region: The latter case corresponds to ω IR = ω UV = Λ and the value of Λ is self-consistently determined by the continuity of the spectral function for a given κ. Thus, this model depends only on κ. In the former case additional considerations are needed to fix ω IR and ω UV , which are described below. We will refer to these two forms as the line model and the step model, respectively.
The NLO result for the spectral function naturally interpolates between the IR and UV regions, but it is not reliable for small ω even at the highest temperature as discussed in section III. However, it can provide some guidance on how to choose ω IR and ω UV . As mentioned above, for ω > T HTL resummation may not be important and the naive and resummed NLO result for the spectral function should agree. As discussed in Appendix B, the resummed and naive NLO results for the spectral function agree well for ω > 2.2 T . Furthermore, the thermal contribution to ρ(ω, T ) is about the same for ω > 2.2 T at the lowest and the highest temperature when normalized by ωT 2 . This indicates that the perturbative calculations are reliable for these values of ω. Therefore, we choose ω UV = 2.2 T . At the highest temperatures, the resummed NLO result is well described by the linear form given by Eq. (19) with κ = κ LO for ω < 0.02 T . Therefore, ω IR = 0.01 T appears to be a reasonable choice. The NLO result for κ is significantly larger than the LO result, implying that the spectral function at low ω is also larger and therefore will match ρ UV (ω, T ) at larger ω. We find that ρ IR (ω, T ) and ρ UV (ω, T ) are equal at around ω = 0.4 T . Therefore, besides ω IR = 0.01 T , we will also use ω IR = 0.4 T and ω IR = 1 T in our analysis.
In Fig. 11 we show the spectral functions obtained from Eq. (20) and Eq. (21) assuming κ = κ NLO in ρ step and ρ line , and three different ω IR at three representative temperatures, T = 1.1 T c , 6 T c and 10 4 T c . From the figure, we see that at the lowest temperature the ρ step (ω, T ) model matches the UV behavior at larger ω without the dip around ω ∼ T of the ρ line (ω, T ) model. The ρ line form with ω IR = 0.01T and ρ step provide an upper and lower bound for the spectral function at T = 1.1 T c . The picture is the same for T = 1.5 T c and T = 3 T c . At T = 6 T c , all forms of the spectral functions provide nearly identical results. At the highest two temperatures, the possible choices of the spectral functions are limited by ρ line with ω IR = 0.01 T and ω IR = T .
Using the models for the spectral functions described above we have calculated the corresponding Euclidean time chromo-electric correlators for different values of κ and compared these with the continuum extrapolated lattice results at each temperature to estimate the heavy quark diffusion coefficient. As discussed in the previous section the continuum extrapolated lattice results need an additional renormalization because the 1-loop renormalization constant, Z E , is not accurate. Therefore, we have matched the correlator obtained from the model spectral function to the continuum extrapolated lattice data at τ T = 0.19. The resulting multiplicative constants C N are slightly different form those shown in Tab. II. This is because the correlators obtained from the model spectral functions are slightly different from G NLO+ E at τ T = 0.19 due to the thermal contribution. We demonstrate this procedure in Appendix B for different model spectral functions. Different forms give different values of κ and this is the dominant source of systematic error in the determination of κ. We have also studied the dependence of κ on the choice of the normalization point in τ and the choice of the renormalization scale. Choosing the normalization point in the range 0.17 ≤ τ T ≤ 0.19 leads to a 8% variation in the resulting κ. Varying the renormalization scale by a factor two results in a similar variation.
Putting everything together we obtain the following estimates for the heavy quark diffusion coefficient from the analysis: although it should be reminded that, as discussed at the end of section III, the lattice data are weakly sensitive to κ at the highest temperature. The dominant uncertainty in the above result comes from the form of the spectral function used in the analysis and the uncertainty of the continuum extrapolated lattice results. We compare our result on κ with the results of other lattice studies [13,[22][23][24][25] in terms of the spatial diffusion coefficient D s which is given by the relation κ/T 3 = 2/(D s T ), in the temperature range T c − 3 T c . This is shown in Fig. 12. We see that our results agree well with the other lattice determinations with the exception of the one in Ref. [13] that is based on charmonium correlators. This is likely due to the fact that the determination of D s from the quarkonium correlators is not accurate since the width of the transport peak is difficult to determine [11,12].
The temperature dependence of the heavy quark diffusion coefficient in the entire temperature region is shown in Fig. 13. We clearly see the temperature dependence of κ 3 /T . The κ obtained on the lattice is not incompatible with the NLO result given the large errors. Inspired by this we fitted the temperature dependence of the lattice result by modeling it on Eq. (16) but keeping the coefficient of m E /T as a free parameter C. From the fit we obtain C = 3.81(1.33), which is larger than the NLO perturbative result C ≈ 2.3302.
We note that our result is significantly larger than the holographic estimate [44]: 2πD s T = 1. Finally, comparing to more experimental quantities, we note that our result for D s at the lowest temperature is much smaller than the calculation from the pion gas [45], which finds, 2πD s T ≈ 17 for T ≈ T c . Experimental determinations of the D-meson azimuthal anisotropy coefficient ν 2 at ALICE [46] and STAR [47] estimate at T ≈ T c κ/T 3 ≈ 1.8 − 8.38 and κ/T 3 ≈ 1.05 − 6.28, respectively. These are in agreement with our findings. All these experimental determinations include mass dependent contributions, while our determination of κ is in the heavy quark limit. Therefore the two should agree up to 1/m corrections.

VI. CONCLUSIONS
In the paper, we have studied the chromo-electric correlator, G E at finite temperature on the lattice with the aim of extracting the heavy quark diffusion coefficient, κ. The calculations have been performed in quenched QCD (SU(3) gauge theory) in order to obtain small statistical errors with the help of the multi-level algorithm. We have studied the dependence of the chromo-electric correlator on the Euclidean time, τ , in a wide temperature range to enable the comparison with weak coupling results. It turned out that the τ -dependence of the electric correlator is poorly captured by the leading order result. Going beyond the leading order result and incorporating the effect of the running coupling in the corresponding spectral function results in a correlation function, G NLO+ E that can capture the τ -dependence of the lattice result much better.
To fully describe the τ -dependence of G E calculated on the lattice, the effect of κ encoded in the low ω part of the chromo-electric spectral function has to be considered. At high ω, we have used forms of the spectral function that are motivated by the next-to-leading order perturbative results. Fitting the lattice results on G E , we have obtained values of κ at different temperatures. We observe that the sensitivity of the chromo-electric correlator to κ is small, varying from few percent at the lowest temperatures to sub-percent at the highest temperatures. This finding is corroborated by a model independent analysis of the chromo-electric correlator, c.f. Figs. 3 and 4. It is this small sensitivity that makes the lattice determination of κ quite challenging. Our main result is summarized in Fig. 13, which shows the temperature dependence of the heavy quark diffusion coefficient. For T < 2 T c our results agree with other lattice determinations, while at higher temperatures they appear consistent with the NLO result. To check to what extent using lattices with aspect ratio N s /N t smaller than four leads to visible finite volume effects we have performed calculations at two temperatures, T = 1.5 T c and T = 10 T c on N 3 s × 12 lattices with N s = 24, 32 and 48. The numerical results are shown in Fig. 14 for some representative values of τ T . As one can see from the figure the finite volume effects are small. We have also attempted to perform an infinite volume extrapolation by fitting the lattice results with a 1/N 3 s form. The corresponding fits are shown in the figure as lines and bands together with the infinite volume result. It is clear from the figure that the differences between the infinite volume result and the lattice results with different N s are of the order of the statistical errors. Therefore, the use of N s = 48 is justified.
As discussed in the main text, to obtain the continuum result for the chromo-electric correlator we first perform the interpolation in τ T and then for each value of τ T we perform the continuum extrapolation using the a/N 2 t form without N t = 12 data or using the a/N 2 t + b/N 4 t form with N t = 12 data included (a and b are fit constants). We have demonstrated this procedure in Fig. 5 for T = 1.1 T c . In Fig. 15 we show this procedure for other temperatures: T = 1.5 T c , 3.0 T c , 6 T c and T = 10 4 T c . We do not show the analysis for 10 T c as it looks similar to the one for T = 10 4 T c . From the figure we see that the slope of the 1/N 2 t dependence increases with decreasing τ T as expected, since the cutoff dependence is larger for smaller τ T . But for the smallest τ T we do not see this tendency. To understand the situation better we show the coefficient of the 1/N 2 t term in the continuum extrapolation as a function of τ T in Fig. 16. We see that the coefficient of the 1/N 2 t term increases in absolute value with decreasing τ T till about τ T = 0.175 and then either flattens off or decreases if τ T is further decreased. We take this as an indication that the continuum limit is not reliable for τ T < 0.175. We also performed continuum extrapolations using lattice data without tree-level improvement and the corresponding results are also shown in Fig. 15 as open symbols. In this case, the continuum limit is always approached from above. The continuum extrapolated result from tree level improved lattice data and the unimproved lattice data agree within errors for τ T ≥ 0.25. In absence of tree-level improvement the continuum extrapolations for smaller τ T are not reliable. In order to understand the main features of the perturbative spectral function corresponding to the chromoelectric correlator at NLO, in Fig. 17 we show the fol-lowing quantity 2(ρ(ω, T ) − ρ(ω) NLO T =0 )/(ωT 2 ) calculated with and without HTL resummation at the lowest and highest temperature. The plotted quantity gives κ in the ω → 0 limit. The naive (unresummed) result is logarithmically divergent at small ω. On the other hand for ω > 2.2 T the resummed and the naive result agree well. This indicates that the NLO calculation is valid in this ω range. We also see that for 2.2 < ω/T < 10 the naive and resummed NLO expressions are negative and their shapes independent of the temperature.
In Fig. 17 we also show the two model spectral functions (line model and step model), where we use κ = 2.75 T 3 for the lowest temperature and κ = 0.088 T 3 for the highest one. At the lowest temperature the step model has a larger finite temperature part than the linear model, while at the highest temperature the opposite is true. The two models also have somewhat different UV behavior. The step model is matched to the zero temperature spectral function and thus ignores the thermal correction in the region 2.2 < ω/T < 10, while the line model incorporates this. The two models thus allow to extract κ using a set of reasonable assumptions about the large ω-behavior of the spectral function.
We match the chromo-electric correlator, G model E obtained from the above model spectral functions at τ T = 0.19 to the continuum extrapolated lattice result to find the optimal value of κ. We demonstrate this procedure in Fig. 18 where we show the continuum lattice result for the lowest and the highest temperature divided by the corresponding G model E . For a given spectral function and the appropriately chosen κ this ratio should be close to one. Since the errors of the continuum extrapolated lattice result are sizable we get a range of κ that is compatible with the lattice result. In Fig. 18 we show the results for κ = 2.75 T 3 and κ = 0.088 T 3 , and for T = 1.1 T c and T = 10 4 T c , respectively. These κ values are chosen to present the extreme values that fit the lattice data within the step model. At T = 1.1 T c , κ = 2.75 T 3 is in the middle of our range for κ (see Eq. (22)), whereas at T = 10 4 T c , κ = 0.088 T 3 is at the upper edge of our range (see Eq. (27)).