Heavy quark diffusion coefficient with gradient flow

We calculate chromo-electric and chromo-magnetic correlators in quenched QCD at $1.5T_c$ and $10^4 T_c$ with the aim to estimate the heavy quark diffusion coefficient at leading order in the inverse heavy quark mass expansion, $\kappa_E$, as well as the coefficient of first mass suppressed correction, $\kappa_B$. We use gradient flow for noise reduction. At $1.5T_c$ we obtain: $1.70 \le \kappa_E/T^3 \le 3.12$ and $1.03<\kappa_B/T^3<2.61$. The latter implies that the mass suppressed effects in the heavy quark diffusion coefficient are 20% for bottom quarks and 34% for charm quark at this temperature.


I. INTRODUCTION
The behavior of a heavy quark moving in a strongly coupled quark gluon plasma (sQGP) can be described by a set of transport coefficients. In particular, the equilibration time of the heavy quarks is related to the heavy quark diffusion. This diffusion can be described as a Brownian motion and, hence, by a Langevin equation that depends on three related transport coefficients [1]: the heavy quark momentum diffusion coefficient κ, the heavy quark diffusion coefficient D s , and the drag coefficient η. In thermal equilibrium these coefficients are related as D s = 2T 2 /κ and η = κ/(2M T ), with M being the heavy quark mass. The heavy quark momentum diffusion coefficient is known in perturbation theory up to mass-dependent contributions at next-to-leadingorder (NLO) accuracy [1][2][3]. We will label this leading term in the T /M expansion as κ E . Moreover, the first mass-dependent contribution when expanding the diffusion coefficient with respect to T /M has been studied in Refs. [4,5], and will be labeled κ B . It is sensitive to chromomagnetic screening and, therefore, is not calculable in perturbation theory [4]. Apart from describing the equilibration time of the heavy quark in a plasma, the diffusion coefficient κ is a crucial parameter entering the evolution equations which describe the out-of-equilibrium dynamics of heavy quarkonium in sQGP [6][7][8].
The NLO correction to κ E is sizable [3], thus calling into question the validity of the perturbative expansion, and inviting instead a strong coupling calculation. Currently, the only analytical strong coupling calculations available are for supersymmetric Yang-Mills theories [9,10], and therefore nonperturbative lattice QCD calculations for the heavy quark diffusion coefficient are heavily desired. However, a direct calculation of the transport coefficients on the lattice can be very challenging, as it involves a reconstruction of the spectral functions from the appropriate Euclidean time correlation functions. The transport coefficient is then defined as the width of the transport peak, a narrow peak at low energy ω. Reconstruction of the spectral function in the presence of a transport peak is a challenging problem, especially since the width of this peak is inversely proportional to the heavy quark mass M [11]. Moreover, the Euclidean time correlators are relatively insensitive to small widths [11][12][13][14][15][16].
The problem of the transport peak can be circumvented by the use of the effective field theory approach. In particular, the heavy quark momentum diffusion coefficient can be related to correlators of field-strength tensor components. The leading contribution in the T /M expansion κ E is related to a correlator of two chromoelectric fields E [6,17], and the T /M correction is related to a correlator of two chromomagnetic fields B [4]. The associated spectral functions ρ E,B (ω) corresponding to these correlators do not have a transport peak, and the heavy quark diffusion coefficients κ E,B are defined as their ω → 0 limit. Moreover, the small-ω behavior is smoothly connected to the UV behavior of the spectral function [17].
The chromoelectric correlator has been calculated on the lattice within this approach in the SU(3) gauge theory in the deconfined phase, i.e., for purely gluonic plasma [18][19][20][21][22], using the multilevel algorithm for noise reduction [23]. It has also been studied out-ofequilibrium using classical, real-time lattice simulations in Refs. [24,25]. During the writing of this paper, the first measurement of the mass-suppressed effects was reported in Ref. [26], also utilizing the multilevel algorithm for noise reduction.
Recently, there has been a lot of interest in using gradient flow [27][28][29] for noise reduction instead of the multi-level algorithm. The gradient flow algorithm is a smearing algorithm that automatically renormalizes any gaugeinvariant observables [30,31] for a sufficiently large level of smearing. The heavy quark diffusion coefficient κ E has been measured with gradient flow in Ref. [32]. Also, preliminary measurements of the chromomagnetic correlator required for κ B have been performed with gradient flow and presented in conference proceedings by two groups [33,34].
In this work, we study both the chromoelectric and chromomagnetic correlators on the lattice using the gradient flow algorithm and determine the diffusion coefficient components κ E and κ B from the respective reconstructed spectral functions. In Sec. II we recall the theory behind the required Euclidean correlators and show the raw lattice measurements of these correlators together with their continuum limits. In Sec. III we then invert the spectral function and provide ranges for κ E,B . The results are summarized in Sec. IV. Preliminary versions of these results have been published in a recent conference proceedings [34].

II. CHROMOELECTRIC AND CHROMOMAGNETIC CORRELATORS
A. Theory background and lattice setup Heavy quark effective theory (HQET) provides a method to calculate the heavy quark diffusion coefficient in the heavy quark limit M πT by relating it to correlators in Euclidean time. The leading-order contribution κ E to the heavy quark momentum diffusion coefficient κ has been expressed in terms of the chromoelectric correlator G E in Refs. [10,17]: where T is the temperature, U (τ 1 , τ 2 ) is a Wilson line in the Euclidean time direction, and E i is the chromoelectric field, which is discretized on the lattice as [17] Recently, the first correction in v 2 to κ, known as κ B , has been put in relation to the chromomagnetic correlator G B [4]: where B i is the chromomagnetic field, which is herein discretized as: For both chromoelectric and chromomagnetic fields, we follow the usual lattice convention and absorb the coupling into the field definition: E i ≡ gE i and B i ≡ gB i . The Euclidean correlators G E and G B are related to the respective heavy quark momentum diffusion coefficient contributions κ E and κ B by first obtaining the spectral functions ρ E,B (ω, T ), where and then taking the zero-frequency limit: The spectral function ρ E (ω, T ) for the chromoelectric correlator G E does not depend on the renormalization in the a → 0 limit, however, the respective spectral function ρ B (ω, T ) for the chromomagnetic correlator G B does [4].
On the other hand, both κ E and κ B are physical observables and the ω → 0 limit of the respective spectral functions does not depend on the renormalization. The two contributions κ E and κ B can then be combined to give the full expression for the heavy quark momentum diffusion coefficient κ [4]: In order to perform the needed lattice calculations, we use the MILC Code [35] to generate a set of pure-gauge SU(3) configurations using the standard Wilson gauge action. The configurations are generated with the heatbath and over-relaxation algorithms, where each lattice configuration is separated by at least 120 sweeps, each consisting of 15-20 over-relaxation steps and 5-15 heatbath steps. We consider two temperatures: a low temperature 1.5T c , and a high temperature 10 4 T c , with T c being the deconfinement phase transition temperature. The temperatures are set by relating them to the lattice coupling β = 6/g 2 0 , which determines the lattice spacing a via the scale setting [36]. This scale setting relates β to a gradient flow parameter t 0 via a renormalizationgroup-inspired fit form, which is then further related to the temperature with T c √ t 0 = 0.2489(14) [36]. For this study, we use lattices with varying numbers of temporal sites, N t = 20, 24, 28, and 34, and with corresponding spatial extents of N s = 48, 48, 56, and 68 sites. Based on our previous study [22], we do not expect there to be a notable dependence on the spatial size of the lattice.
To measure the Euclidean correlators we rely on the gradient flow algorithm [27][28][29]. The Yang-Mills gradient These equations are an explicit representation of Adapting these equations for a pure-gauge lattice theory with link variables gives us the differential equatioṅ where V τF are the flowed link variables. We choose S Gauge to be the Symanzik action. The lattice simulations with N t = 20, 24, and 28 are evaluated numerically with a fixed step-size integration scheme [29], while the N t = 34 lattice is evaluated with an adaptive step-size implementation [37,38]. For further analysis, we need the data points from all lattices at the same flow time positions. Therefore, we use cubic spline interpolations with simple natural boundary conditions in order to provide the data along a common flow-time axis. The full list of parameters and statistics are given in Table I. The gradient flow evolves the unflowed gauge fields A µ to the flowed fields B µ , which have been smeared with a flow radius √ 8τ F . This smearing systematically cools off the UV physics and automatically renormalizes the gauge-invariant observables [31]. This renormalization property of the gradient flow is especially useful for the correlators G E,B , which otherwise require renormalization on the lattice. The renormalization of the correlators can be calculated in lattice perturbation theory, like in Ref. [39], but the lattice perturbation theory may have poor convergence [40]. In previous multilevel studies of the chromoelectric correlator [21,22], a perturbative oneloop result for the chromoelectric field renormalization Z E was used [39]. As gradient flow automatically renormalizes gauge-invariant observables, such a factor Z E is not needed in this study, as has been observed already in the previous studies of G E using gradient flow [32]. The continuum-and flow-time-extrapolated result for G E at 1.5T c obtained using gradient flow agrees with the continuum extrapolated results obtained using the multilevel algorithm and with the one-loop result for Z E at the level of a few percent, indicating that for the β range considered in the calculations of G E , the perturbative renormalization is fairly accurate. Moreover, in a recent lattice study of a different but similar operator, where a chromoelectric field was inserted into a Wilson loop [41], it was shown explicitly that Z E → 1 at sufficiently large flow times. For chromomagnetic fields, renormalization is required both on the lattice and in continuum [5]. As the renormalization property of gradient flow is generic to all gauge-invariant observables [31], the chromomagnetic correlator should require no additional renormalization on the lattice either.
On the other hand, since the gradient flow introduces a new length scale √ 8τ F , we have to make sure it does not contaminate the measurements at the length scale of interest τ -the separation between the field-strength tensor components. The most basic condition for ensuring that the flow has enough time to smooth the UV regime, while preserving the physics at the scale τ , would be a < ∼ √ 8τ F < ∼ τ /2. The upper limit of this condition was further restricted in Ref. [32], by inspecting the LO perturbative behavior of the flow [42], to be (τ − a)/3. In our experience, slightly larger flow times are still fine; hence, we use a slightly relaxed limit: Moreover, we note that instead of dealing with the scales √ 8τ F and τ separately, the relevant scale for these Euclidean correlators is in fact the ratio of the scales, √ 8τ F /τ . This can be inferred from the leadingorder result of the chromoelectric correlator at finite flow time [42]: where ξ 2 n = x 2 n T 2 /τ F and x n = τ + n/T . Likewise, the early NLO result from Ref. [43] also shows an affinity to this ratio. Using the units of √ 8τ F /τ , the condition of suitable flow times from Eq. (14) becomes We use these limits for both G E and G B . In order to reduce the discretization errors further, we define a tree-level improvement by matching the LO continuum perturbation theory result [17], to the LO lattice perturbation theory result [19], We then define a tree-level improvement of G E as [32] G imp where the improvement is restricted to zero-flow-time discretization effects, because the lattice perturbation theory result for Symanzik flow is not known. We use the same tree-level improvement for the chromomagnetic correlator as for the chromoelectric correlator, since in the continuum limit these correlators are the same at leadingorder. We also used the clover discretization for the chromomagnetic correlators in addition to the one given in Eq. (4). The clover discretization was used in Ref. [26]. We check that in the continuum limit, the clover discretization and the one given by Eq. (4) yield identical results within errors. The corresponding analysis is discussed in Appendix A. This fact gives us confidence that the discretization errors are well under control.

B. Lattice measurements
In Fig. 1, we present both electric and magnetic correlators of the raw lattice data, normalized with Eq. (17) and tree-level improvement at different flow times for a single representative lattice size N t = 28. We observe the statistical errors decreasing as the ratio √ 8τ F /τ increases, and that for √ 8τ F /τ > 0.1 the curves at different flow times seem to converge toward a common shape. This shape seems to be shared between both G E and G B .
Next, we perform the continuum extrapolations of both correlators. First, we interpolate the data for each lattice in τ T at a fixed flow-time ratio with cubic spline interpolations. Since the correlators G E,B are symmetric around the point τ T = 0.5, we set the first derivative of the splines equal to zero at τ T = 0.5. We perform a linear extrapolation in 1/N 2 t = (aT ) 2 of the correlators at the fixed interpolated τ T , and fixed flow-time ratio positions, using lattices N t = 20, 24, 28, and 34 for large separations τ T > 0. 25. For small separations τ T < 0.25, we drop the N t = 20 lattice from the extrapolation. As an example, we show the continuum extrapolations at different τ T and √ 8τ F /τ in Fig. 2. The χ 2 /df of the continuum extrapolation is around 1 or smaller. For small √ 8τF/τ for the Nt = 28 lattice at T = 1.5Tc. We see that with increasing ratio, the correlators converge toward a common shape across the whole τ T range.
τ some continuum extrapolations have large χ 2 , indicating that the cutoff effects are too large to obtain reliable results. We also perform continuum extrapolations including a 1/N 4 t term for lattices with N t = 16, which corresponds to a O(a 4 ) continuum extrapolation. These continuum extrapolations agree with the ones shown in Fig. 2 within errors. Further details on the continuum extrapolations are discussed in Appendix A.
We present the continuum limits at the edges of the √ 8τ F /τ range, within witch we will later take the zeroflow-time limit, and see that the continuum values vary less when √ 8τ F /τ is changed than when τ T is changed. Hence, the thermal effects of heavy quark diffusion dominate the shape of these correlators. Figures 3 and 4 show the final continuum limits of G E and G B , respectively, for both measured temperatures as function of τ T . Similarly to what we observed in Fig. 1, both correlators exhibit similar behavior at fixed temperatures according to the shape and order of magnitude. As mentioned above for the chromomagnetic correlator, we also perform calculations using clover discretizations and verify that the same continuum limit is obtained for this discretization. This is shown in Appendix A.   To further inspect this similarity, in Fig. 5 we plot the ratio G E /G B along the fixed √ 8τ F /τ axis, and observe a near-constant behavior toward large separations τ T . From here, we can already deduce that the contribution to the heavy quark diffusion coefficient from the chromomagnetic correlator G B is only going to differ from the contribution of the chromoelectric correlator G E by less than 5%. In Fig. 3, we also show the zero-flow-time limit for the chromoelectric correlator G E , which will be discussed further in the next section.

A. Modeling of the spectral function
We now turn to extracting κ E and κ B from G E and G B , respectively, using Eqs. (5)    for modeling the spectral function closely follows the approach laid out in our previous work [22]. This approach uses the perturbative information on the spectral function at large ω, where this information is expected to be reliable. For both correlators, the spectral function ρ E,B is known at the NLO level [26,44]. We chose to model the spectral function such that in the UV regime at zero flow time it follows the T = 0 part of the NLO spectral function; however, we chose the scale so that the NLO part vanishes, leaving us with only the LO part [17]: The coupling has been evaluated at the five-loop 1 level in the MS scheme at the scale µ opt ω , which for ρ E reads [44],  For the electric spectral function ρ E , we further change to the NLO EQCD scale [45] when ω ≈ T or smaller. As we discuss below, the LO or NLO result for the UV part of the spectral function is not accurate, and we hence multiply it by a normalization factor C n to take into account higher-order corrections, i.e., we perform the replacement ρ LO E,B → C n ρ LO E,B . A similar normalization constant was used in the analysis of Refs. [21,32]. The determination of C n is discussed at the end of this subsection. For ρ E , we do not model the flow-time dependence of the spectral function, as one is able to take the zero-flow-time limit before the spectral function inversion.
For the magnetic spectral function ρ B , the situation is more complicated due to the required renormalization [4,5]. In order to study the chromomagnetic correlator renormalized correlator in the MS scheme: where h 0 is a constant and γ 0 = 3/(8π 2 ) is the anomalous dimension of the chromomagnetic field [26]. In principle, the renormalization constant Z flow can be calculated in perturbation theory; however, in practice we know from our previous calculation [22] that the NLO perturbative results are not reliable enough to fully describe the lattice data. Hence, Z flow is fixed by comparing the perturbative result to the lattice result on G B . Using the NLO result from Ref. [26] and neglecting the distortions due to finite flow time (i.e., setting h 0 to zero), Eq. (25) gives a flowtime-dependent UV part of the chromomagnetic spectral density: where β 0 = 11/(16π 2 ) is the leading coefficient of the β function, and As with ρ E , we choose the scale µ opt in such a way that Eqs. (26) and (22) are equal up to a constant, Z flow : As in the case of the chromoelectric spectral function, here we make the replacement ρ UV to take into account higher-order perturbative corrections. The determination of the normalization constant C n is discussed below, and now C n will also contain the unknown normalization factor Z flow .
The perturbative spectral functions described so far cover the UV regime of our model spectral functions. Alone, these UV spectral functions would give κ E,B = 0, and hence an infrared contribution needs to be added leading to finite κ E,B . We note that while in general the chromomagnetic spectral function depends on the renormalization scheme (MS, gradient flow, etc.) and scale, its low-frequency limit does not since κ B is a physical quantity. This has been shown explicitly in weak-coupling calculations [4]. One can work with the physical (RGinvariant) chromomagnetic spectral function by scaling out the anomalous dimension, or one can equally well work with the chromomagnetic correlation function in the gradient flow scheme at some finite, but sufficiently small, τ F , and extract κ B .
In order to extract the κ E,B , we then follow the procedure laid out in our preceding study [22] and model the spectral function with a family of Ansätze. For the large-ω regime in the UV, we assume the LO perturbative spectral function at T = 0 as ρ UV from Eq. (22) to hold. while for small ω in the IR, the spectral function is given by We assume that ρ E,B (ω, T ) = ρ IR (ω, T ) for ω < ω IR and ρ E,B (ω, T ) = ρ UV E,B (ω, T ) for ω > ω UV , where ω IR and ω UV are the limiting values of ω for which we can trust the above behaviors. In the region ω IR < ω < ω UV , the form of the spectral function is generally not known, and this lack of knowledge will generate an uncertainty in the determination of κ E,B . Hence, for a given value of κ E,B , we construct the model spectral function that is given by ρ UV E,B in ω > ω UV , ρ IR E,B in ω < ω IR , and a variety of forms of ρ E,B (ω) for the intermediate ω IR ≤ ω ≤ ω UV , such that the total spectral function is continuous. For the functional forms of the spectral function in the intermediate ω values, we consider two possible forms based on simple interpolations between the IR and UV regimes: where θ(ω) is a step function. The case described in Eq. (31) corresponds to ω IR = ω UV = Λ with the value of Λ self-consistently determined, i.e., the value of Λ is set by requiring the model spectral function to be continuous. We will refer to these two forms as the line model and the step model, respectively. In our previous analysis, we determined that the NLO spectral function takes the linear form for ω < 0.02T , and converges to the UV form at ω > 2.2T , and hence we use the same ω IR = 0.01T and ω UV = 2.2T as in Ref. [22] for the line model (30) and for both chromomagnetic and chromoelectric spectral functions. The correlation functions obtained from the model spectral functions through Eq. (5) will be labeled as G model E,B . The spectral representation of G E,B given by Eq. (5) also holds at finite lattice spacing (a = 0) and finite flow time (τ F = 0), as long as the spectral function ρ E,B is replaced by a lattice equivalent ρ lat E,B (a, τ F ). The spectral function ρ lat E,B (a, τ F ) only has support for ω < ω max . In the case of meson correlators, a similar ρ lat has been explicitly constructed in the free case [46]. In this work, point-like meson sources and sinks are used. One often uses correlation functions of extended meson operators to improve the signal of the ground state. The spectral function of such extended meson correlators has also been calculated in the free theory [47]. It was found that for extended meson operators, the support of the spectral function shifts toward smaller ω values [47], and their shape is modified at large ω but not at small ω [47]. The operators obtained from gradient flow can be viewed as extended operators, and therefore the shape of the corresponding spectral functions at large ω will be different compared to the unsmeared case, and the support of the spectral function will shift toward smaller ω. However, the small-ω limit of ρ lat E,B (a, τ F ) will not depend on a or τ F to a good approximation, because the correlator G E,B is not sensitive to a or τ F , provided that τ a and τ √ 8τ F . Therefore, in principle, one can extract κ E,B even at finite a and τ F . However, as this is valid only for ω < ω max , the large-ω part of ρ lat E,B (a, τ F ) cannot be described by the continuum perturbative result. We model the UV part of the spectral function ρ UV E,B with Eq. (22), up to a multiplicative constant. The difference between these and the continuum spectral functions is not expected to be large in terms of the correlators G E,B at τ T > 0.25, which is the relevant τ range for the determination of κ E,B . The above statement about the dependence of the spectral function at large ω on the flow time appears to contradict the perturbative analysis of Ref. [32]. However, we note that in Ref. [32] the analytic continuation was done in terms of the Matsubara frequency, while here we consider continuation in terms of τ : t → −iτ . These two methods of analytic continuations lead to different results, unless the spectral function decays like 1/ω 2 for large ω. Finally, we note that the cutoff effects in ρ lat are not limited to the large-ω region. There are cutoff effects proportional to (aT ) 2 ∼ 1/N 2 t which affect ρ lat at all values of ω. However, these are quite small for the N t > 16 used in our calculations.
So far, we have presented the correlators G E,B normalized with Eq. (17), which assumes a constant coupling. In Fig. 6 we include the running coupling in the analysis, and divide the continuum limit of the correlators for √ 8τ F /τ = 0.239 with Eq. (5), using Eq. (22) with the scales (23) and (28) for ρ E and ρ B , respectively. The corresponding correlators are labeled as G LO+ E,B . We see from Fig. 6 that with this normalization the τ dependence of the corresponding ratios is greatly reduced. In particular, at the highest temperature 10 4 T c , only very little τ dependence can be seen for τ T ≥ 0.25. The τ depen-dence observed for τ T < 0.25 is most likely due to the fact that our continuum extrapolation is not reliable at such small τ [22]. Thus, a large part of the τ dependence of G E,B comes from the running of the coupling constant. On the other hand, the values of the ratios G E,B /G LO+ E,B differ significantly from one, even at relatively small τ . A similar trend for G E /G LO+ E was observed in Ref. [22]. It was speculated in Ref. [22] that the fact that G E /G LO+ E is roughly a constant that is different from one may be due to the one-loop renormalization of the lattice correlator not being reliable. However, as discussed in Sec. II A, the one-loop renormalization of the chromoelectric correlator is quite reliable. This leads us to the conclusion that the NLO results for the spectral function may not be reliable and an additional normalization constant, C n , has to be introduced as an extra fit parameter. The normalization constants C n are shown in Appendix B. For the chromoelectric correlator, the normalization constant C n is very close to the one obtained in our study using the multilevel algorithm [22]. In the case of the chromomagnetic correlator, the constant C n also contains the unknown matching between the gradient flow scheme and the MS scheme, as mentioned before. We suspect that the fact that the NLO result can describe the lattice correlators at small τ only up to a constant C n is due to the presence of the Wilson line and the Polyakov loop in the definition of the correlators. These do not contribute at order g 4 , but will start contributing at higher orders. It is also known that the weak-coupling result for the Polyakov loop only works at temperatures T > 5 GeV [48]. At higher orders, the presence of the Wilson line and the Polyakov loop most likely changes the overall normalization of the correlator, but not its τ dependence.

B. Flow time dependence of the correlators
To get rid of distortions due to gradient flow, the lattice results for G E should be extrapolated to zero flow time. The limit to zero flow time has to be taken after the continuum limit to avoid the large discretization effects at small flow times. Also, it was argued in Refs. [33,43] that the inversion of the spectral function via Eq. (5) is mathematically well defined only in the zero-flow-time limit. As discussed in the previous subsection, it is possible to generalize the spectral representation in Eq. (5), for nonzero lattice spacing and flow time, if the corresponding spectral function only has support for ω values smaller than some ω max .
We also note that in lattice studies of shear viscosity, spectral function inversion at finite flow time has given satisfactory results [49,50]. Moreover, in recent studies of latent heat, it has been observed that the order of the continuum and zero-flow-time limits can be switched as long as one is careful to only take the limits in regimes where the functional forms used are justified [51]. Therefore, we present our main analysis following the conventional order continuum limit → zero-flow-time limit → 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 spectral function inversion for the main analysis, but we also present an analysis where these steps are taken in a different order. To perform the extrapolation to the zero-flow-time limit, we use a linear ansatz in τ F . A linear behavior is expected, as the small-τ F behavior is just a leading correction to the τ behavior due to flow. Moreover, for the chromoelectric correlator G E , the linear behavior has been seen at the NLO level of perturbation theory [43]. Starting with the chromoelectric correlator G E , we present examples of linear zero-flow-time extrapolations at a few chosen τ T values in Fig. 7. As expected, we see a clear linear dependence in the range where the extrapolation can be performed. We observe that the correlator G E decreases with increasing flow time. The whole range of τ T dependence of the zero-flow-time results was already presented in Fig. 3. As one can see from that figure, the flow-time dependence is not very large in the considered flow-time window. In particular, the shape of the correlator does not change significantly with the flow time and it is very similar to the shape of the correlator extrapolated to zero flow time. Thus, the determination of κ E is not significantly affected by 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 the nonzero flow time. Therefore, one can also model the spectral function corresponding to nonzero flow time and determine κ E . The effects of the small residual distortion of the correlator due to gradient flow on κ E can be taken care of by performing a zero-flow-time extrapolation for the resulting κ E . This analysis strategy will be discussed in the next subsection.
The flow-time dependence of the chromomagnetic correlator is shown in Fig. 8, and appears to be quite different from the flow-time dependence of the chromoelectric correlator. The flow time dependence of G B appears to be roughly linear, but its slope has the opposite sign. This difference is expected and probably comes from the nontrivial renormalization of G B , cf. Eq. (25). This renormalization is taken care of at leading-order in G LO+ B . Normalizing the chromomagnetic correlator by G LO+ B , instead of by G norm , largely reduces the flow-time dependence. This is shown in Fig. 9. In the case of the chromomagnetic correlator we do not take the zero-flowtime limit, but instead model the spectral function for nonzero flow time using Eqs. (26), (29), (30), and (31), and then perform the zero-flow-time extrapolation of κ B obtained from this modeling, as will be described below.

C. Results: κE
The extraction of κ E proceeds as follows. We take the continuum and zero-flow-time limit data and perform a least-squares fit to Eq. (5) with either of the models of the spectral function ρ E,B (ω). In addition to having κ E,B as a fit parameter, we also enforce a normalization in the fit by finding a normalization coefficient C n (τ min T ) fit parameter such that G E (τ min T )/G model E (τ min T ) = 1. To estimate the contributions from the systematic errors, we perform these fits with different values of τ T min , vary the scale µ of the running coupling by a factor of 2, and perform the fit with either the line (30) or step (31) models of the spectral function. To get the final estimate for the heavy quark momentum diffusion coefficient, we then take the full spread of the subset of these fits for which the ratio G E (τ T > τ min T )/G model E (τ T > τ min T ) = 1 is within 1.5σ.
For the G E data, which was first extrapolated to the continuum limit and then to the zero-flow-time limit, this procedure gives for κ E at T = 1.5T c and at T = 10 4 T c . The T = 1.5T c result gives a slightly improved range for κ E compared to our previous multilevel study [22], which had 1.31 < κ E /T 3 < 3.64. Although slightly smaller, it is also in agreement with the other existing results for this temperature: 2.31 < κ E /T 3 < 3.70 from Ref. [33], 1.8 < κ E /T 3 < 3.4 from Ref. [21], 1.55 < κ E /T 3 < 3.95 from Ref. [20], and 1.3 < κ E /T 3 < 2.8 from Ref. [26]. The T = 10 4 T c result is in agreement with our previous result 0 < κ E /T 3 < 0.1 [22]. The new result has slightly larger errors due to the gradient flow analysis having more strict fit regimes; however, we can for the first time observe a nonzero minimum for κ E /T 3 at very large temperature. Both of these κ E values can be reexpressed as a position-space momentum diffusion coefficient D s = 2T 2 /κ [17] as: 0.64 < D s T < 1.17 for T = 1.5T c and 12.5 < D s T < 100 for T = 10 4 T c .
We now turn to a question of how the result for κ E depends on the order of the limits. First, in Fig. 10 we show the extracted values of κ E /T 3 both at the zero-flowtime limit and at a finite flow time for both the line (30) (filled symbols) and the step (31) (empty symbols) models for the spectral function ρ(ω). Only points that are within the regime 0.2 < √ 8τ F /τ < 0.3 where reasonable zero-flow-time extrapolation can be performed and that satisfy the condition G E (τ T > τ T min )/G model E (τ T > τ T min ) = 1 within 1.5σ are shown. In addition, we show in different colors the different choices of the normal- ization point τ T min . We observe that the variation between the models is the dominant source of error and that the variation within the flow time is small in comparison. Moreover, with reasonably high τ T min ≥ 0.275, we have enough data points to perform a linear zero-flowtime extrapolation, which we can see agrees closely with the results we get from data that has been extrapolated to zero flow time before spectral function inversion, although with much larger errors. If we were to do the κ E extraction purely at finite flow time, the full variance due to the different fit forms would give us 1.5 ≤ κE T 3 ≤ 3.2 for T = 1.5T c and 0.007 ≤ κ E /T 3 ≤ 0.18 for T = 10 4 T c . Therefore, the variance for a given finite flow time is much larger than the difference between the continuum extrapolated κ E extractions.
Furthermore, we inspect whether it matters that the continuum limit is taken before everything else, as has been done so far. If we were to instead extract the κ E at finite lattice spacing and then take the continuum limit as a linear extrapolation of the extracted κ E values, we would get the result in Fig. 11. We see that the continuum limit of the κ E extracted at finite lattice spacing replicates both the zero-flow-time result, and the results at finite ratio √ 8τ F /τ . Hence, all results presented above would remain unchanged even if the continuum limit had been taken last, because of the large uncertainties in κ E due to the modeling of the spectral function.

D. Results: κB
We now turn to the chromomagnetic correlator G B and extraction of the respective κ B . Based on the above analysis for κ E , we can safely assume that one can get a very good estimate of the zero-flow-time-extrapolated value even when limiting the analysis to a finite flow time. Our analysis strategy here closely follows the case of the 12. Spectral function of chromomagnetic correlators obtained from the two fit forms described in the text at two different representative flow-time ratios. Only the mean value is shown and the statistical errors are hidden for better visibility.
chromoelectric correlator. First, we fix the normalization constant C n and then vary κ B to obtain the best agreement of the lattice correlator with the model correlator.
To demonstrate this point for the step Ansatz, we write where Λ T ∼ T is some IR cutoff. We treat κ B as a fit parameter, while C n (τ F ) is adjusted such that G model B from the above equation exactly matches the continuum lattice result for G B at τ = τ min . The values of C n are shown in Appendix B as a function of τ F . In Fig. 12 we show the spectral function corresponding to the chromomagnetic correlators for √ 8τ F /τ = 0.25 and 0.3. As one can see from the figure, the flow-time dependence of the spectral function is rather mild.
The flow time behavior of the extracted κ B /T 3 is shown in Fig. 13, where again the filled symbols show the extraction using the line ansatz (30), the empty symbols show the extraction with the step model (31), and different colors depict the different choices of τ T min . We observe less curvature in the extracted κ B values than we saw for κ E in Fig. 10. If we take the total variation at finite flow time to be the error of κ B , we get for T = 1.5T c that 1.23 < κ B /T 3 < 2.74. We can then proceed to take the zero-flow-time limit in the linear regime √ 8τ F /τ ≥ 0.25, similar to what we learned to work with in the case of G E . In the zero-flow-time limit, we get the final result for κ B :

IV. CONCLUSIONS
In this paper, we have studied the chromoelectric and chromomagnetic correlators in quenched QCD with the aim to determine the heavy quark diffusion coefficient, including the subleading correction in the inverse quark mass. We used gradient flow for noise reduction and showed how to control the distortions due to nonzero flow time in the calculations of the transport coefficients κ E and κ B . To obtain the heavy quark diffusion coefficient, we used a parametrization of the spectral functions that relies on the NLO result at large energies, and smoothly matched to the expected linear behavior at small energies. The effects of the nonzero flow time can be incorporated into the high energy part of the spectral function. We verified this in the calculations of κ E , where we obtained κ E from the chromoelectric correlator extrapolated to zero flow time, as well as by calculating an effective κ E from the chromoelectric correlator at finite flow time and then extrapolating to zero flow time. Our main results are summarized in Eqs. (32), (33) and (35). Our results for κ E agree with the previous determinations [22,33,52] within the estimated uncertainties. The value of κ B we obtained agrees with the very recent result obtained using the multilevel algorithm and nonperturbative renormalization based on the Schrödinger functional [26]. We have seen that the dominant uncertainty in the determination of κ E and κ B comes from the modeling of the spectral functions at low energies. Using the lattice results for v 2 from Ref. [13] for charm and bottom quarks v 2 charm 0.51 and v 2 bottom 0.3 (c.f. Fig. 6 of Ref. [13] where v 2 th = v 2 /3 was shown), we estimate that the mass-suppressed effect on the heavy quark diffusion coefficient is 34% and 20% for charm and bottom quarks, respectively.
The extraction of the heavy quark diffusion constant strongly relies on using the NLO result for the spectral function at large energies. It is assumed that the NLO result can describe the τ dependence of the correlators up to a multiplicative constant. To test this assertion further, it would be desirable to perform calculations at larger N τ , so that reliable continuum extrapolations are possible for smaller values of τ . Another way to obtain more reliable continuum-extrapolated results is to use the Symanzik-improved gauge action. We plan to implement such an improved analysis in the near future. Finally, once the full one-loop perturbative matching between the MS scheme and the gradient flow scheme at small flow times becomes available, we will redo our analysis by converting to the MS scheme and taking the zero-flow-time limit. In this appendix, we discuss discretization effects and continuum extrapolations in more detail. In Fig. 14 we show different continuum extrapolations for the chromoelectric correlators as a function of τ for a few representative values of the flow time. We perform extrapolations assuming a 1/N 2 t form for the discretization errors and vary the range in N t , and also include a 1/N 4 t term in the continuum extrapolations with N t = 16 lattices. For τ T > 0.25 different continuum extrapolations agree well with each other.
The χ 2 /df of the continuum extrapolations are shown in Fig. 15. For τ T > 0.25, different continuum extrapolations of the chromoelectric correlators agree well, and the χ 2 /df of the continuum extrapolation is close to one or smaller. For smaller τ T the χ 2 /df is large, indicating that the continuum extrapolations are not reliable. We perform a similar analysis for the chromomagnetic correlators. Some results are shown in Fig. 16  As discussed in the main text for the chromomagnetic correlators, we use two discretization schemes: the simplest one given by Eq. (4), which can be labeled as the corner discretization, and the clover discretization which was also used in Ref. [26] (c.f. Eqs. (2.2)-(2.4) herein). These two discretizations must agree in the continuum, but could lead to quite different result at nonzero lattice spacings. As a result, the tree-level improvements for these two discretization schemes are also different. The leading-order result for the clover discretization without the g 2 and C F factors has the form whereq andq are given by Eqs. (19) and (20) respectively. We use this to implement the tree level improvement for the clover discretization scheme. In Fig. 17 we show the continuum limit of the flowed chromomagnetic correlator with the clover discretization and the tree-level improvement (A1). We show results for two different flow times, for which the expected 1/N 2 t behavior can be clearly seen in the lattice data in both cases. We also compare the continuum extrapolated results obtained with the corner and clover discretizations, and the corresponding tree-level improvements, in Fig. 18. As one can see from the figures, the continuum results obtained with the two discretization schemes are in excellent agreement. The tree-level improvement reduces the discretization effects and therefore, aids robust continuum extrapolations. However, as discussed in Ref. [22] it is not necessary if the lattice spacing is sufficiently small or, equivalently, if N t is large enough. Small lattice spacings are needed for reliable continuum extrapolation at small τ T . If τ T is not very small, the continuum extrapolation can be performed without tree-level improvement [22]. To check to what extent our conclusions on the continuum result of the chromomagnetic correlator depend on the tree-level improvement, we perform continuum extrapolations of the chromomagnetic correlator with the corner discretization scheme but using the "wrong" tree-level improvement, namely, the tree-level improvement for the clover discretization. The corresponding continuum results are also shown in Fig. 18  labeled as "corner+cloverNorm". For τ T < 0.35 we see small, but statistically significant, differences compared to the continuum results obtained with proper tree-level improvement, but for larger values of τ T the tree-level improvement is not essential for reliable continuum extrapolations.

Appendix B: Normalization parameter
For completeness, we also show in Figs. 19 and 20 the normalization coefficient C n for both G E and G B respec-tively. We observe that C n has a very mild dependence on the flow time. This can be used as an indication that modeling ρ lat E,B with the running coupling version of the leading-order ρ E,B is reasonably well motivated. The C n values for the chromoelectric correlator are well in agreement with the ones we reported in our preceding study [22]: ∼ 1.73 for T = 1.5T c and ∼ 1.2 for T = 10 4 T c . The C n for the chromomagnetic field is slightly larger than the respective factor for G E .