Heavy quark momentum diffusion from the lattice using gradient flow

We apply the gradient flow on a color-electric two-point function that encodes the heavy quark momentum diffusion coefficient. The simulations are done on fine isotropic lattices in the quenched approximation at $1.5\,T_c$. The continuum extrapolation is performed at fixed flow time followed by a second extrapolation to zero flow time. Perturbative calculations of this correlation function under Wilson flow are used to enhance the extrapolations of the non-perturbative lattice correlator. The final estimate for the continuum correlator at zero flow time largely agrees with one obtained from a previous study using the multi-level algorithm. We perform a spectral reconstruction based on perturbative model fits to estimate the heavy quark momentum diffusion coefficient. The approach we present here yields high-precision data for the correlator and is also applicable for actions with dynamical fermions.


I. INTRODUCTION
The Yang-Mills gradient flow has proven to be a powerful tool for gauge theories since it was proposed in [1][2][3][4][5]. Its applications on the lattice can be useful in several ways, for example reducing high frequency background noise from gauge configurations or setting the physical scale in lattice QCD [6][7][8]. Observables calculated under gradient flow are claimed to be automatically renormalized at sufficiently large flow time if the continuum limit is taken [3,9]. This is a convenient property when dealing with operators whose renormalizations are troublesome, for instance the energy momentum tensor, which has already been studied under flow [10][11][12][13]. Even more applications of the gradient flow can be found in [14].
Studies on transport coefficients like the heavy quark momentum diffusion coefficient, whose knowledge is of phenomenological interest, also benefit from the gradient flow. Non-perturbative determinations of such coefficients require high precision estimates of corresponding two-point correlation functions whose signals are often overshadowed by high frequency background noise even at large Monte Carlo sample sizes. Techniques like the multi-level algorithm [15] and link-integration [16] can ameliorate the problem but are not applicable for actions with dynamical quarks. For now the gradient flow seems to be the only way to obtain high precision data of correlation functions in simulations with dynamical quarks. Discussions on how to flow dynamical quarks and their applications in dynamical QCD can be found in [5,13,17].
In this paper we want to study how the well-known correlation function of color-electric fields which contains information about heavy quark transport behaves under Yang-Mills gradient flow and how to perform its continuum and flow-time-to-zero extrapolation. This colorelectric correlator [18] has been intensively studied in perturbation theory [18][19][20][21] and also non-perturbatively [22][23][24][25]. For this study we restrict ourselves to pure SU (3) gauge theory at a temperature of T ≈ 1.5 T c . We compare our results with a previous non-perturbative measurement of the correlator from [24] that was obtained in the same setting but with the multi-level algorithm as a noise reduction method. This allows us to cross-check the results obtained from the gradient flow approach. The extrapolated correlator in the continuum limit at zero flow time is then used for spectral function reconstruction, for which we use multiple theoretically motivated fitting models, similar to what was done in [24]. Finally we compare the heavy quark momentum diffusion coefficient obtained from the extracted spectral function with those from other studies [24,26].
The paper is structured as follows: we start by recalling the definition of the gradient flow and how it is realized on the lattice (section II). Afterwards we briefly present some perturbative calculations of the color-electric correlator at non-zero flow time and how we can use them to enhance the non-perturbative lattice data (Section III). Section IV is then devoted to the analysis of this data and how we carry out the double extrapolation. The final estimate for the renormalized continuum correlator at zero flow time is then used to obtain the heavy quark momentum diffusion coefficient through spectral reconstruction in section V. We summarize our findings in section VI.

II. GRADIENT FLOW
The Yang-Mills gradient flow evolves the gauge field to some flow time τ F along the gradient of the gauge action. The flowed gauge field B µ (τ F , x) is defined through

arXiv:2009.13553v1 [hep-lat] 28 Sep 2020
where A µ is the ordinary gauge field. In order to fulfill gauge invariance one can easily construct the flowed covariant derivative and field strength following the standard Yang-Mills gauge theory: In perturbation theory, a composite local operator O(x, τ F ) consisting of gauge fields can be expanded in τ F as a superposition of the renormalized operator defined at zero flow time [4]: Here c i are some coefficients that can be calculated perturbatively. In order to obtain the renormalized operator in the original theory one needs to invert Eq. 3 and go back to zero flow time. We will show how to carry out this procedure in practice in section IV D.
To leading order a perturbative calculation in D dimensions [3] provides a simple explanation for the effect of the flow on the gauge field: The flow averages the gauge field with a local Gaussian kernel. The mean-square radius of the Gaussian distribution is √ 8τ F and is sometimes called the "flow radius." On the lattice we need to use the discretized counterparts of the continuum flow equations. Specifically, we use a Symanzik improved version of the flow called Zeuthen flow [27], which eliminates O(a 2 ) cutoff effects that are present in the ordinary Wilson flow. We solve Eq. 1 numerically using an adaptive step-size Runge-Kutta method for Lie groups [28].
How much flow can be applied to a gauge field before measuring a correlation function on the lattice? Consider the correlator of two operators whose temporal separation is τ . On the one hand, the flow radius should be large enough to suppress lattice discretization effects, and on the other hand, we need to make sure that the separation between operators is larger than the flow radius, so that the correlation between them is not contaminated. This leads to a general allowed flow time range of We can make the latter limit more rigorous by performing a perturbative determination of the flow radius where a specific correlator starts to suffer contamination due to the overlap of the operator smearing radii. We will do this for the color-electric correlator below.

III. THE COLOR-ELECTRIC CORRELATOR IN PERTURBATION THEORY
The correlator that we want to study under flow was first proposed in [18] and is related to the thermalization rate of a heavy quark in a hot plasma. It arises from the infinite quark mass limit of the meson current-current correlator. The advantage of studying this specific correlator is that the spectral function encoded in it has no sharp transport peak, in contrast to the perturbative behavior for current-current and stress-stress correlators. As a result, the relevant transport coefficient is much easier to obtain. The correlator is a product of electric fields that sit on a Polyakov loop, which is why it is referred to as the color-electric correlator. It is defined as [18] where β = 1 T is the temporal extent which is equal to the inverse temperature, and U (τ 2 , τ 1 ) is a Wilson line in the Euclidean time direction. E i is the color-electric field in the geometrical normalization, which we discretize following [18]: Here, U i (τ, ⃗ x) is an SU(3) link variable in the i th direction.
The leading-order perturbative calculation in continuum of the correlator under gradient flow was performed in [20]. We perform the analog calculation on the lattice here in order to illustrate the behavior of the correlator under Wilson flow compared to its continuum counterpart.
The LO perturbative lattice correlator under Wilson flow reads where I 0 , I 1 are modified Bessel functions of the first kind and C F = (N 2 c − 1) 2N c = 4 3 is the quadratic Casimir of the fundamental representation. Details of the derivation of this result can be found in the appendix A 1. We divide the perturbative correlator by g 2 C F and equip it with the superscript "norm" since we will later use it to normalize our non-perturbative lattice measurements.
In Figure 1 we show an example for different flow times in the continuum and on a lattice with N τ = 20. The vertical dashed line in this figure is the lower limit of the distance for which the flowed continuum correlator deviates by less than 1% from its non-flowed counterpart and is given by τ ≳ 3 √ 8τ F [20]. For the lattice version of this limit one needs to shift the allowed distance by about one lattice spacing, and together with Eq. 5 one obtains a serviceable flow interval of Because of this limitation one can only hope to extract the large-distance behavior via the gradient flow method. Fortunately this is also the region that carries the information about the heavy quark momentum diffusion coefficient. From Figure 1 we can see that the correlator grows as τ −4 at small separations, and therefore spans several orders of magnitude in value. In order to improve visibility we can remove this dominant behavior from the non-perturbative lattice data that we present later by normalizing it with the perturbative leading-order correlator. The normalized data is also convenient for the extrapolations as its values are all of the same order of magnitude.
Correlators calculated on the lattice suffer from cutoff effects. One can remove the leading-order contribution of these effects from the non-perturbative lattice data by utilizing the perturbative lattice and continuum correlators (see Eq. 8 and [20]). This technique is known as tree-level improvement [29]. In appendix A 2 we explain in detail how to perform the improvement for the color-electric correlator under flow. Note that in our study the non-perturbative correlator is obtained using Zeuthen flow, while the perturbative lattice calculations are done with Wilson flow. The main point of appendix A 2 is that it is reasonable to improve the Zeuthen-flowed correlator by using the non-flowed leading-order lattice correlator. Since we want to normalize the data with the leading-order continuum correlator, we end up with the dimensionless ratio (cf. Eq. A6 and A7) on which the continuum and flow-time-to-zero extrapolations will be performed.

IV. FLOWED COLOR-ELECTRIC CORRELATOR ON THE LATTICE
In this section we present our measurements and analysis of the color-electric correlator on fine isotropic lattices at T ≈ 1.5 T c . We perform the continuum extrapolation at fixed flow time and subsequently the extrapolation to τ F = 0.

A. Lattice setup
A summary of the parameters of the gauge configurations used in this study can be found in Table I. All configurations are generated in the quenched approximation using heatbath and overrelaxation updates. The gauge action is the standard Wilson action. After thermalization (5000 heatbath sweeps), we save a configuration after every 500 combined sweeps, where one such sweep consists of one heatbath and four overrelaxation sweeps. The lattice boundary conditions are periodic for all directions.
In order to compare our correlator measurements with previous results from [24], we use the same lattice dimensions and β-values to set the temperature to about 1.5 T c . The lattice spacing a is determined via the Sommer scale r 0 [29] with parameters taken from [30] and updated coefficients from [31]. We use the state-of-theart value r 0 T c = 0.7457(45) [30].
For statistical mean and error estimation of observables measured on the configurations we perform a bootstrap analysis with 10000 samples for each lattice.

B. Cutoff effects and signal-to-noise ratio under flow
Calculations of the color-electric correlator without gradient flow exhibit huge statistical errors even for the large amounts of gauge configurations used in this study (see Table I). Beyond the first few lattice spacings the obtained data is almost pure noise at zero flow time. By increasing the flow time, high-frequency fluctuations of the gauge fields are suppressed and the signal-to-noise ratio of the correlator is expected to increase. Additionally, cutoff effects are expected to decrease with increasing flow time.
The minimum amount of flow that is needed to produce renormalized operators is √ 8τ F T ≈ aT = N −1 τ . Via Eq. 9 one can convert this minimum flow time into a lower bound for the separation τ T . Because of this we henceforth exclude the N τ = 16 lattice (see Table I) from the analysis. The limits are then dictated by the N τ = 20 lat-  [29] with parameters taken from [30] and updated coefficients from [31]. We use r0Tc = 0.7457(45) [30]. The coarsest lattice (Nτ = 16) does not enter in the continuum extrapolation.
tice, and so we have to flow up to at least √ 8τ F T = 0.05, and we can only obtain renormalized correlator data for τ T ≥ 0.2. These limits are also applied to the finer lattices and will carry over to the continuum and flow-time-tozero extrapolations.
In Figure 2 we show the tree-level improved correlators on all lattices at one intermediate ( √ 8τ F T = 0.06) and one larger flow time ( √ 8τ F T = 0.10). The dashed vertical line depicts the lower boundary for the distance resulting from the upper boundary for the flow time (Eq. 9). As expected, the flow drastically improves the signal and ameliorates cutoff effects, but it also contaminates the correlation at smaller distances if too much flow is applied. The continuum values shown in these figures will be explained in the following section. Figure 3 provides a more detailed look into how the correlator behaves under flow at fixed separations on the finest available lattice (N τ = 36). After the initial "signal acquisition phase" there is only a minor effect that the flow exerts on the longer-distance correlation. The second lattice separation is particularly interesting since it is large enough to prevent immediate operator contamination but small enough to inherit almost no noise at zero flow time. At this separation we can see an initial rising behavior under flow, which we attribute to the renormalization caused by the flow. (This behavior can be observed on all lattices of Table I; the rising behavior is not visible for the higher separations as here the dominant effect of the flow is the improvement of the signal.) A similar renormalization-induced behavior can also be observed for the electric correlator in perturbative lattice QED under Wilson flow, for which we explicitly calculate the "tadpole-type" next-to-leading order contributions in appendix B. Our understanding of the behavior in The grey error bands indicate the underlying data; to reduce clutter they are hidden if statistical errors would be larger than 5%. In retrospect we only consider the flow times for the flow extrapolation where the data points are shown explicitly (see section IV D). The dashed vertical line is the minimum amount of flow that is needed for the coarsest considered lattice (Nτ = 20, see Table I). The horizontal positions of the various markers depict the flow limits from Eq. 9 (with a = 1 (Nτ T ), Nτ = 20).
plete and the correlation then behaves almost linearly for some time, before it is then substantially contaminated. After a continuum extrapolation we will exploit the seemingly linear region to extrapolate the correlator to zero flow time. The dominant flow-time-dependent contributions in this region will essentially consist of gradually accruing contact terms of the operator. In principle these exist as soon as the flow time is non-zero, and their contribution is enhanced by the fact that the operator is not truly local due to the discretization of the electric fields (Eq. 7).

C. Continuum extrapolation
The physical correlation function is at vanishing lattice spacing and flow time, so we need to perform a double extrapolation. But it is important to perform the extrapolations in the right order: the continuum extrapolation should be performed first, and the flow-time-to-zero extrapolation should be carried out on the continuum correlation functions. Otherwise, lattice-spacing issues which are ameliorated by flow re-emerge in the flow-time-tozero extrapolation.
The continuum extrapolations should only be performed at the separations of the finest lattice that is available, as the information is most reliable there. All data from coarser lattices need to be interpolated to these separations, for which we use cubic splines. For a given lattice and flow time we only interpolate inside the interval given by Eq. 9, with the exception of the one data point that lies just outside of the lower separation limit. By also including this slightly contaminated point we can always interpolate at the lower limit without explicitly using the contaminated data, which extends the usable flow range. At the left boundary of the spline the second derivative is set to zero (natural condition) and at the middle point (τ T = 0.5, which in practice is the right boundary) the first derivative is set to zero (symmetry). We perform the interpolations on every bootstrap sample and take the bootstrap mean and standard deviation of the interpolated points for statistical error estimation.
Once the interpolated values are obtained, the continuum extrapolation is carried out, independently for each distance τ and flow time τ F , by performing a weighted fit on the data with the Ansatz with fit parameters m and b = G cont τ,τ F G norm cont τ,τ F =0 . The Wilson action that we use to generate the gauge configurations has leading discretization errors of order a 2 , which explains this choice. We estimate the statistical error by performing the weighted fit on substitute bootstrap samples that we draw from a Gaussian around the bootstrap mean of the interpolations with their error as the width. The bootstrap mean and standard deviation of the substitute samples then serve as the statistical estimates. In

D. Flow-time-to-zero extrapolation
The continuum correlators at fixed flow time still need to be extrapolated back to τ F = 0. Because of the limited knowledge we have of the functional form of the flow effect we use the simple linear Ansatz where m and b = G cont τ G norm cont τ are tau-dependent fit parameters. For statistical error estimation we proceed as in the continuum extrapolation.
Based on the continuum correlators we estimate by hand for which separations and flow times this Ansatz is reasonable (see Figure 5 and also Figure 3): For the lower boundary of the flow time the criterion is the statistical precision of the data. There are two ways to improve the precision: 1. apply more flow, 2. increase the number of gauge configurations. For the larger distances where one needs more configurations for the same statistical accuracy (in comparison to the smaller distances), the minimum amount of flow from Eq. 5 is not enough to give sufficiently small statistical errors on the data. Therefore we only use continuum data points in the flow-time-to-zero extrapolation that are flowed enough so that they have a maximum relative error of 1% ⋅ τ T . This value is chosen to yield a lower boundary for the flow time where a linear Ansatz is reasonable. This also implies that the minimal separation for which we can obtain renormalized correlator data is shifted. For the upper boundary of the flow time we are limited by Eq. 9 and by the number of data points of the coarsest lattice; if there are less than two non-contaminated data points for any given lattice, we cannot perform the interpolation  Table I) in comparison with revised continuum correlator from multi-level method (perturbatively renormalized to NLO [32]) using lattice data from [24]. and continuum extrapolation. For the coarsest lattice this turns out to be at a flow time of √ 8τ F T > 0.133.
With all of these limitations we can reliably extrapolate the renormalized continuum correlator in the range τ T ∈ [8 36, 0.5]. In Figure 6 we show the final renormalized continuum correlator at zero flow time. The errors shown here only depict statistical uncertainties. We also carry out a new continuum extrapolation of the lattice correlators from a previous study [24] which used the multi-level algorithm and link integration techniques to obtain a signal for the color-electric correlator. The old method only works in pure gauge theory and needs to be renormalized approximately with the help of perturbation theory [32]. Note that the statistics used in [24] are much smaller than the ones in this study (Table I), especially for the finer lattices.
The two results agree in their overall shape, while the slight offset can be explained by missing higher-order contributions to the renormalization of the multi-level result, the systematic uncertainty introduced by the flow extrapolation, and the difference in lattices and statistics.
A note on systematic uncertainties: In principle, the correlator data shown in Figure 6 is subject to systematic uncertainties through scale setting, auto-correlation between gauge configurations and cross-correlation between the correlator distances, the non-flowed tree-level improvement of the Zeuthen-flowed correlators, and finally the interpolation and extrapolation Ansätze.

V. HEAVY QUARK MOMENTUM DIFFUSION COEFFICIENT A. Spectral function of the color-electric correlator
It was shown in [18] that the heavy quark momentum diffusion coefficient can be obtained from where ρ(ω) is the spectral function of the color-electric correlator, which is related to the Euclidean correlator we compute here through the integral relation In the following we will calculate the leading-order continuum spectral function under gradient flow, which will lead to a discussion about the Kramers-Kronig relation and flowed retarded correlators. The (continuum, zero flow time) leading-order spectral function is well known [21]. But we did not find a study of the spectral function for the Euclidean function at finite flow depth, so we will consider this problem here. We start from the leading-order correlation function at flow time τ F , where Γ(s, x) = ∫ ∞ x dt t s−1 e −t is the incomplete Gamma function and d the number of space dimensions. Fourier transforming the correlator to frequency space using the (periodic) Kronecker delta function ∫ β 0 dτ e iknτ = βδ kn,0 leads tõ Setting the dimension to its physical value d = 3 and performing the analytic continuation we find the spectral function to be The real parts of the incomplete Gamma functions do not vary with their incomplete argument. The value of the real part of the incomplete Gamma function is therefore the same as the ordinary Gamma function. The leading order spectral function is which is the the same result as the non-flowed spectral function calculated in [21]. Therefore, at leading perturbative order, the spectral function for the flowed and for the non-flowed Euclidean correlation function are the same. At first this appears remarkable; two distinct Euclidean functions reconstruct to the same spectral function, which implies that Eq. 14 must not apply to the flowed correlation function. The connection between the retarded correlator and the spectral function is made by the Kramers-Kronig relation (see [33]). The real and imaginary parts of a meromorphic function, which is analytic in the closed upper half-plane, are related via Kramers-Kronig integral formulas if the function vanishes fast enough for a growing complex argument. The mathematical requirements of the function can be translated to the necessity of causality of the underlying physical theory.
The introduction of gradient flow, that is, replacing local operators with smeared-out relatives, breaks energymomentum conservation through new contact terms. The commutator correlator of smeared-out operators no longer vanishes outside the light cone, thus breaking causality. Alternatively, we can appeal to Eq. 4, where e −τ F p 2 represents an exponential suppression for Euclidean 4-momentum p, but is an exponential enhancement for timelike Minkowski 4-momenta. Mathematically speaking, the exponential suppression along the Euclidean time axis τ will turn to exponential growth after Wick rotation. This happens during the analytic continuation of the Euclidean to the retarded correlator, which no longer fulfills the requirements needed for the Kramers-Kronig relation, thereby breaking its connection to the spectral function.
This discussion proves that spectral reconstruction can only be attempted on data which has been extrapolated to zero flow time (as well as to the continuum). As we have already emphasized, this extrapolation should be performed after extrapolating to the continuum limit. With the extrapolated data we can apply spectral reconstruction techniques to estimate the heavy quark momentum diffusion coefficient κ, which we will do in the next section.

B. Extraction of the heavy quark momentum diffusion coefficient
There are many approaches for a spectral reconstruction [34][35][36][37][38][39]. Here we will proceed similarly to [24] and reconstruct the spectral function based on theoretically motivated model fits of the correlator data obtained from the gradient flow method that is shown in Figure 6. The fit ansatz is a combination of two parts: one is valid in the small frequency regime (ω ≪ T ) and the other is valid in the large frequency regime (ω ≫ T ). In the intermediate regime an interpolation is required. The small frequency part can be described by infrared asymptotics [18], and the large frequency part can be calculated in perturbation theory with the help of asymptotic freedom [21,40], or, when considering the higher order correction, C F , r 20 , r 21 , can be found in [21,24] and g 2 , a s are evaluated using four-loop running [41]. Eq. 21 and Eq. 22 will be treated equally in our fits as explained in [24]. To account for the corrections when incorporating φ IR and φ UV , which are valid only in their own regimes, we need to introduce two trigonometric functions, e (α) n (y) ≡ sin(πny), e (β) n (y) ≡ sin(πy) sin(πny), (23) where y ≡ x 1+x and x ≡ ln 1 + ω πT . Putting all this together yields three types of models as in [24]; however, in this work we only adopt the most theoretically justified one, where µ ∈ {α, β}, i ∈ {a, b}. Substituting the true spectral function in Eq. 14 with the model spectral function, we obtain a model correlator G model with which we define a weighted sum of squares for a fit: Here δG cont (τ ) is the statistical error of the lattice data G cont (τ ) which we show in Figure 6. Now it is clear that the parameters to fit are κ T 3 and c n . As in [24], two strategies are considered. The first strategy is just a normal fit and in the second strategy we constrain that at the maximum frequency the spectral function reproduces its UV asymptotics. In strategy I we do not stop after 200 iterations as in [24] as we have seen that this criterion makes almost no difference. For both strategies in our fit routine we set the absolute tolerance for χ 2 to 0.0001 and the maximum number of iterations to 2000. In our fits we limit ourselves to n max ∈ {4, 5} and we use all the data points shown in Figure 6, that is, τ T ≥ 8 36.
In Figure 7 we show an example of the fit correlators and spectral functions when using model ρ correlators are also similar. This also holds for all the other choices of models or fit strategies and we obtained The resultant κ T 3 is summarized in Figure 8. For each point the error bar is the statistical uncertainty we obtain from a bootstrap analysis and we take the spreading over various models as the systematic uncertainty of our estimation. Finally, based on the central values of Figure 8 we obtain a range for κ T 3 :  [26]. With the estimated range for κ T 3 we can also estimate the heavy quark diffusion coefficient in the heavy quark mass limit (M ≫ πT ) according to D = 2T 2 κ [18]:

VI. CONCLUSION
We calculated the color-electric correlator on four isotropic lattices at T ≈ 1.5 T c in the quenched approximation under gradient flow. We found that the gradient flow technique significantly improves the signal of the correlator on the lattice. It is necessary to perform both a continuum and a flow-time-to-zero extrapolation; we show that the continuum limit must be taken first, and we perform both extrapolations on our data within the range τ T ∈ [8 36, 0.5]. We found that an Ansatz linear in 1 N 2 τ and one linear in τ F is suitable for the respective extrapolations. The overall shape of the correlator obtained via the flow method in the a → 0 and τ F → 0 limit largely agrees with its multi-level counterpart in the a → 0 limit. And just like the shape of the correlator, the heavy quark diffusion coefficient obtained from the flowed correlator is also consistent with but slightly larger than that from the correlator computed with the multi-level approach, according to the χ 2 fits performed on the correlators using theoretically well-established models.
This work affirms the validity of the gradient flow technique when performing the continuum and flow-time-tozero extrapolations and provides a method for extending these studies to full QCD.
where ξ and α 0 are the gauge and flow gauge fixing parameters and we used the tilde short-hand notatioñ In the following we will use both Feynman and flow Feynman gauge. The leading-order expression reads which is analogous to the continuum expression but exchangingk → k. The spatial integration can be performed analytically where I 0 and I 1 are modified Bessel functions of the first kind. Inserting this into Eq. A2 and introducing a Schwinger parameter in order to rewrite the denominator leads to the expression given in Eq. 8. contribution of these effects by comparing the perturbative continuum and lattice correlator ("tree-level improvement" [29]). This can be achieved by defining the improved distances τ T through which are used to shift the correlator values according to Another possibility is given by multiplying the correlator data at their original distances with ratio of the leading-order continuum and lattice correlator: In principle it should not matter whether one improves the distances or the correlator values themselves. However, the interpolations of the correlator data that are necessary for the continuum extrapolation will change slightly. The data will not be distributed evenly if one improves the separations, and one obtains slightly smaller large separations (e.g. τ T = 0.5 → τ T ≈ 0.48) which is unwished-for since most of the information about the lowfrequency part of the spectral function is contained at the largest τ T . Consequently, we choose the straightforward strategy of Eq. A5 to perform the tree-level improvement.
In the end we want to normalize our lattice measurements to the leading-order continuum correlator and so we end up with In principle it is necessary to take both correlators on the right side of Eq. A6 at the same flow time τ F ; however, currently we have only calculated the leading-order color-electric correlator under Wilson flow (Eq. A2), which is substantially different from the more complicated Zeuthen flow [27] that we use for the numerical simulations.
In Figure 9 we compare the Zeuthen-flowed correlator improved by the non-flowed leading-order correlators (left panel) with the Wilson-flowed correlator improved by the Wilson-flowed leading-order correlators at the same flow time (right panel) on a coarse lattice with N τ = 16. While the behavior under flow of these two quantities differs greatly, their naive flow-time-to-zero extrapolations using a linear Ansatz (see also section IV D) differ only by about 1%, which is shown in Figure 10. This difference is mainly caused by the discretization, and we expect it to be even smaller for finer lattices, which is why we decide to use the non-flowed tree-level improvement on the Zeuthen-flowed measurements: The continuum and flow-time-to-zero extrapolations are performed on this improved correlator normalized by the leading-order continuum correlator. for the second smallest distance and beyond if the statistical errors are small enough. (The rising behavior is not visible for the higher separations as here the dominant effect of the flow is the improvement of the signal.) Figure 11 shows the Wilson-flowed (not Zeuthen) colorelectric correlator as a function of flow time on a N τ = 16 lattice where the effect is clearly visible. We expect this behavior to stem from the renormalization of the colorelectric field which is caused by the flow. The purpose of this section is to give a qualitative explanation for this behavior. This is done by calculating the two-point function of the electric field in compact QED lattice perturbation theory in vacuum as a function of the flow time to next-to-leading order: This enables us to qualitatively analyze the renormalization behavior of the electric field under gradient flow with a smaller number of diagrams compared to QCD. Nonetheless, the renormalization properties should be the same in both theories. The leading order contribution has already been derived in the thermal theory (Eq. A2) and so we find At next-to-leading order in the lattice spacing a there are three diagrams, which arise from the expansion of the electric field operator (Figure 12), the action ( Figure 13) and the Wilson flow equation (Figure 14).
The tadpole diagram ( Figure 12) is a lattice artifact and describes the mixing between the electric field operator E and higher-order operators such as E 3 which are induced by the lattice (lattice tadpole corrections). The contribution of this unphysical type of diagram vanishes in the continuum limit as there is no equivalent diagram in the continuum. However, for non-continuumextrapolated calculations the contribution of the diagram needs to be considered and the interpretation is already well understood and we refer to [42] for a detailed analysis. One of the electric fields is expanded up to O(a 4 ) leading to 12. Next-to-leading order diagram of the expansion of the color-electric field operator. On the left the diagram is drawn in our notation; below the horizontal dotted line are the non-flowed propagators and vertices, and the vertical axis above the dotted line represents evolution in flow time. On the right side the same diagram is drawn in the notation introduced in [4]. We will denote this diagram to be contribution (a). This diagram is a pure multiplicative correction to the leading order and describes the tadpole renormalization of the operators. The integral can be calculated analytically yielding where I 0 denotes the modified Bessel function of the first kind. The term −3 2 a 4 I 2 (τ F ) is shown in Figure 15. From the figure it is clear that for zero flow time the diagram has a significant correction to the leadingorder behavior of the correlator. However, if more flow is applied, this contribution vanishes. One can interpret this behavior as follows: the lattice introduces highdimension operators with coefficients containing powers of a: E → E + a 4 E 3 . The flow smears out the fields and introduces a smooth momentum-cutoff, so that the expectation value of such high-dimension operators is suppressed by an inverse power of the flow time; the effects a 4 E 3 are parametrically of order a 4 τ F 2 . Another NLO effect arises from interaction terms in the expansion of the action. The lowest order interaction is a four-point vertex. The corresponding diagram is the first diagram contributing to the photon self-energy up to this order in a and is shown in Figure 13. It has the same meaning as in the continuum: the renormalization of the coupling. Since the loop momentum lives at zero flow time, the corresponding correction is expected to be independent of the flow time. Consequently, the correlator as a function of gradient flow takes the form where the integral I 1 is defined as This shift is exactly the leading-order contribution to the difference between the bare and renormalized coupling. The last diagram ( Figure 14) emerges from the nonlinear part of the flow equation leading to a four-point vertex at an intermediate flow time τ ′ . It is a pure gradient flow correction to the photon self-energy. In QED this is a pure lattice artifact; in QCD similar terms also appear in the continuum theory but the lattice introduces a distinct gauge-invariant subclass of effects with no continuum analogue. The unique NLO correction in lattice QED is where L β (Y, τ ′ ) is the Lagrange multiplier field needed to introduce the Wilson flow equation in its Lagrangian form, see [4] for more details. We note that in this case the integral of the leading-order contribution is not recovered, instead an additional power ofK 2 emerges. Since this diagram represents a pure Wilson flow effect, it is not surprising that I 3 (τ F ) vanishes for zero flow time.
The integral as a function of the flow time is shown in Figure 15. From the figure we can learn that the contribution of this diagram asymptotically reaches a constant value for large flow times. This is essentially a shift in the flow time that one can absorb into the definition of the flow time. To see this, we consider the case that more flow loops are attached to the propagator connecting the electric fields. These additional loops will lead to a correction of the form LO(τ, τ F ) 1 + c 1 (τ F )a 2K 2 (B7) the sort shown in Figure 15 lead to c n (τ F ) = (c 1 (τ F )) n , and the corrections resum into the form e +K 2 f (τ F ) , which can be interpreted as a finite order-a 2 renormalization of the flow time. This effect has been observed before in previous lattice studies in [43,44].
Combining the contributions of the diagrams, we can now write down the full correlation function at next-toleading order: Evaluating this expression leads to an initial rise of the QED correlator as a function of flow time due to the renormalization of the electric field. We consider this effect to be present in the QCD color-electric correlator ( Figure 11) as well, originating from similar renormalization effects as found in the NLO calculation in QED.