Lattice $B$-field correlators for heavy quarks

We analyze the color-magnetic (or"$B$") field two-point function that encodes the finite-mass correction to the heavy quark momentum diffusion coefficient. The simulations are done on fine isotropic lattices in the quenched approximation at $1.5\,T_c$, using a range of gradient flow times for noise suppression and operator renormalization. The continuum extrapolation is performed at fixed flow time followed by a second extrapolation to zero flow time. Perturbative calculations to next-to-leading order of this correlation function, matching gradient-flowed correlators to MS-bar, are used to resolve nontrivial renormalization issues. We perform a spectral reconstruction based on perturbative model fits to estimate the coefficient $\kappa_B$ of the finite-mass correction to the heavy quark momentum diffusion coefficient. The approach we present here yields high-precision data for the correlator with all renormalization issues incorporated at next-to-leading order, and is also applicable for actions with dynamical fermions.


I. INTRODUCTION
Heavy-ion collisions create a new state of matter, the quark-gluon plasma [1][2][3].Studying this state of matter is complicated by the fact that most participants-quarks and gluons-undergo multiple interactions throughout the collision, ending in hadronization.Therefore, their final yields and momentum distributions encode the early dynamics in a complicated and indirect way.Hard probes are objects with enough energy and/or momentum to provide more direct information about the early stages of the collision.Principle among these are heavy quarks, whose evolution in the quark-gluon plasma has been extensively studied both experimentally [4][5][6] and theoretically [6][7][8].
Most heavy quarks are produced at relatively low transverse momentum, meaning γv ≲ 1.In this regime, they interact with the medium through a series of scatterings which can be collectively described as momentum diffusion and drag [7,9,10].Close to rest, the momentum-diffusion coefficient κ is determined by a force-force correlation function which can be rewritten as a color-E-field color-E-field correlation function in the medium, with the adjoint color-electric field group-theory factors connected by a fundamental Wilson line [11].The associated spectral function is the analytical continuation of the thermal Euclidean correlation function of two E fields along a fundamental Polyakov loop, normalized by the trace of the Polyakov loop [12].This opens the possibility of computing this quantity on the lattice, something several groups have pursued [13][14][15][16][17].
Because the charm quark is not extremely heavy, and because some charm and bottom quarks are created with somewhat larger momenta, it may also be important to consider finite-velocity corrections to this momentumdiffusion picture.Bouttefeux and Laine showed [18] that the next order in velocity is determined by a colormagnetic-magnetic correlator on a Wilson line, which is the continuation of a magnetic-magnetic correlator along a fundamental Wilson line in the thermal Euclidean path integral.This quantity can also be investigated on the lattice, and several groups have already tried to do so [19][20][21].However, treating this problem is more complicated, as Laine already showed [22]: the BB correlator on a Wilson line renormalizes in a nontrivial way, which must be taken into account both in fitting lattice data and in interpreting the result in terms of heavy-quark diffusion.
The goal of this paper is to investigate the correlation function of two magnetic fields on a Wilson line with the use of lattice QCD and gradient flow.We hope to improve on previous studies in three ways.First, we perform a high-statistics study using a wide range of very fine lattices in order to control continuum and small-flow extrapolations.Second, we treat correlations within the data analysis process more carefully than previous studies (including our previous study of electric field correlations).For instance, results at a fixed separation and nearly equal flow depths are highly correlated, whereas many analyses have treated different flow depths as if they are independent; similarly, the small-flow extrapolation for nearly equal separations is also highly correlated.Third, we treat the issues of operator renormalization more carefully when carrying out the finite flow time extrapolation.Here we take advantage of a recent next-to-leading order (NLO) result for the effect of gradient flow on magnetic field correlators along a Wilson line to provide a fully NLO matching between gradientflowed lattice results, perturbative calculations in the modified minimal subtraction scheme (MS scheme), and the physical diffusion coefficient.The use of these relations greatly improves the extrapolation over gradient flow depth, eliminating a logarithmic dependence which would otherwise arise.
An outline of the paper is as follows.In Sec.II we present the theoretical framework of this study, including the definition of lattice object, analytical continuation, gradient flow and a matching procedure for converting lattice results at finite gradient flow depth to the physical values.Next, Sec.III describes the details of our lattice setup and the evaluation of magnetic-field correlation functions.Sec.IV provides the details of the matching procedure.Sec.V describes the continuum and zero-flow limits.This section represents a main innovation in our treatment, since the renormalization of the operators as a function of gradient flow depth is fully incorporated in our extrapolations, and since data correlations across flow times and between separations is taken into account.Then Sec.VI presents the analytical continuation of our results to a Minkowski spectral function.Finally, Sec.VII presents our results and conclusions.

II. RENORMALIZATION OF MAGNETIC-FIELD CORRELATORS UNDER GRADIENT FLOW
The basic object we shall consider is the correlation function of two magnetic field operators along a Wilson line.In Euclidean space this will be a fundamentalrepresentation Polyakov loop, and the value should be normalized by the mean trace of the Polyakov loop: Here β = 1/T is the inverse temperature, U (t 1 , t 2 ) is a fundamental-representation Wilson line along the time direction, G 12 is the color magnetic field strength under the geometrical normalization convention, and L is the trace of the Polyakov line.All operators occur at a common spatial point, only time is varying.The analytical continuation of this object to real time is the Wightmanordered correlation function of two magnetic field operators with a Wilson line running up the real time axis, between the operators, down the real-time axis, and across the Euclidean branch: The associated spectral function ρ(ω) is related to G B (τ ) in the standard way, In a completely analogous way one can also define G E , the correlation function for electric fields, where the − sign accounts for two factors of i in the continuation of Minkowski to Euclidean electrical fields. 1he two spectral functions ρ E and ρ B then determine the momentum-diffusion coefficient κ via where ⟨v2 ⟩ is the mean-squared velocity of the heavy quark.The velocity, as determined by thermodynamics, depends on the ratio T /M , which is the way the dependence on the thermal heavy quark mass M enters.The operators in Eq. ( 1), Eq. ( 4) require regularization and their regularized values will generically depend on the scale or procedure used.For the case of electric field correlators on a Wilson line this regularization turns out to be harmless; the operator determining G E turns out to be finite and regulator-independent within MS and gradient flow regulators, and its value in these two schemes is equal.But the issue is more complicated for G B (τ ), as we review below.
The lattice itself provides a regularization which could be used.But the regularization calculation has not been carried out to our knowledge, and in general latticeregularized operators tend to possess large renormalizations and poorly convergent perturbative expansions, so this approach is not attractive.In any case, achieving a high signal to noise from a lattice simulation will generally require the use of gradient flow.Gradient flow automatically provides renormalized operators with finite correlation functions, and it eliminates any detailed dependence on the lattice spacing or lattice implementation.But it is important to correctly understand what this regularization actually does in terms of the flow time dependence of G B .
Briefly, gradient flow is a procedure for transforming gauge-field backgrounds A µ (x) within the path integral, removing their short-distance fluctuations.One defines a flowed variable B µ (x, τ F ) through a boundary condition and an evolution equation, where the superscript B indicates that the covariant derivative and field strength are those constructed using the flowed field B µ (x, τ F ), not the original field A µ (x).
For more discussion see Refs.[23][24][25][26].After applying this procedure to a flow depth τ F , the gauge fields are free of short-distance fluctuations beyond approximately the momentum scale and all correlation functions are rendered finite.We will write the magnetic field correlation function evaluated on the flowed fields as G flow,µF It is the quantity which one directly measures on the lattice.The quantity G E has a finite τ F → 0 limit which is approached with O(τ F /τ 2 ) corrections.However, as we will see, G flow,µF B also contains logarithmic-in-τ F corrections which must be correctly handled.
Bouttefeux and Laine [18], followed by Laine [22] considered G B in the MS renormalization scheme and its relation to the physical diffusion coefficient for heavy quarks.They found the relation between the physical diffusion coefficient and the MS-regulated and renormalized 3 correlation function G MS,μ B to be (11) Here g 2 γ 0 = g 2 C A /8π 2 = 3g 2 /8π 2 is the leading-order anomalous dimension for a magnetic field operator inserted along a temporal Wilson line (with C A = 3 the adjoint Casimir) with the gauge coupling g 2 = 4πα s .We assume that g 2 = g 2 MS (μ).Since G B involves two B-field insertions, the correction has an additional overall factor of 2 relative to that in Ref. [22].This expression implies that G MS,μ B depends explicitly on the MS renormalization scale μ through a Callan-Symanzik equation: It remains to relate G MS,μ

B
with what we actually measure on the lattice, which is G flow,µF B .Since the latter is independent of the μ scale, its relation with G MS,μ B must also involve logarithms of μ, balanced by the gradient flow scale µ F .Very recently, de la Cruz et al. have found the relation [27] between these quantities, valid to leading order in small τ F /τ 2 , This result was also recently derived by Brambilla and Wang [28].The μ dependence in Eq. ( 13) and in Eq. ( 11) cancel, rendering a combined result which is independent of μ at the NLO level.Besides these renormalization effects, we expect that G flow,µF B will contain additional polynomially suppressed 3 Meaning that 1/ε factors have been subtracted off.
flow effects, that is, effects of order τ F /τ 2 .To extrapolate these away, we will be interested in small values of τ F , satisfying τ F T 2 ≪ 1 and therefore µ F ≫ T .If we use a single MS scale μ in applying both Eqs. ( 11) and ( 13), this will lead to large logarithms in at least one equation.We should then expect α 2 s ln 2 (large) effects at the next order.Such uncontrolled higher-order logarithms can be avoided, as usual, by evaluating Eqs. ( 11) and ( 13) at two different MS scales, each chosen appropriately to avoid large logarithms in the matching, and evolving G MS,μ B between these scales using the Callan-Symanzik equation Eq. (12).This also requires an accurate determination of the scale-dependent gauge coupling.Combining these ideas, the following approach appears appropriate to obtain the renormalized color-magnetic correlator G phys.B : 1. Determine the MS gauge coupling g 2 MS (μ) = 4πα s at one scale, μref , accurately, using the approach of Refs.[25,29] (more details in Sec.IV).From this, use the standard MS beta function to determine the gauge coupling over the desired range of scales.
2. For each relevant4 flow time, measure G flow,µF B (τ F ) on the lattice and take its continuum limit at fixed flow time.

Match the continuum-extrapolated
with the help of Eq. ( 13).This requires choosing a matching scale μτ F for μ.This could be, for example, the typical scale for gradient flow effects, or the scale where G flow,µF In any case we end up with G MS,μτ F B . 4. To avoid large logarithms in Eq. ( 11), run G MS,μτ F B to the "physical scale" μT using Eq. ( 12), incorporating the running of the coupling using a two-loop or higher-loop beta function.Two sensible choices for μT would be μT,LO ≡ 2πT, the typical scale for thermal physics, and the scale where In any case we end up with G MS,μ T B .
In the event that either Eq. ( 11) or Eq. ( 13) involves a relatively large correction, we can also attempt to incorporate higher-loop multiple-log contributions by the substitution (1 + x) → exp(x), where x = γ 0 g 2 . . . is the NLO correction.A summary of steps 3-5 is then where we have introduced the total matching factor Z match .The resulting quantity, G phys.

B
(τ F ), must then be extrapolated to zero flow time, assuming that the flow time dependency is of form τ F /τ 2 respectively.
In the following we will employ this procedure in our data analysis to correctly account for the renormalization of G B , such that the spectral reconstruction can be performed using G phys.

B
. The approach involves a somewhat arbitrary choice for the scale μT and for the scale μτ F .The dependence on these choices formally cancels at the NLO level, but it will introduce residual nextto-next-to-leading order (NNLO) effects.By considering two somewhat different choices for each scale, we will determine the residual renormalization-point dependence of our results, which we will include in our systematic error budget.

III. LATTICE DETAILS
The color-magnetic field correlators are measured on four isotropic lattices at temperature T ≃ 1.5T c , where T c is the confinement/deconfinement phase transition temperature determined via the Sommer parameter r 0 [30] in Ref. [31]. 5 The gauge configurations are generated with the standard Wilson gauge action using the heat bath and overrelexation algorithms.To ensure our calculation is performed on fully thermalized configurations, the first 4000 sweeps, each consisting of one heat bath update and four overrelexation updates, have been discarded.conditions are used for all directions.The lattice spacing a is determined via the r 0 scale [31,32].The finite temperature lattice setup is summarized in Table I.For the calculation of the coupling constant, we use additional zero-temperature lattice simulations at the three finest lattice spacings considered here, see Table II.Before measuring the color-magnetic field correlator on the lattice, we first evolve the gauge fields to the desired flow time range using the Symanzik-improved gradient flow (Zeuthen flow) [33].
The generation of gauge configurations and the correlator measurements, as well as the gradient flow evolution is performed using the SIMULATeQCD suite [34][35][36].The correlator data analysis and spectral reconstruction is carried out using our software toolkit correlators flow [37].

IV. DETERMINATION OF GAUGE COUPLING AND MATCHING FACTOR
In order to make use of Eq. ( 18) and the matching procedure explained in Sec.II, the gauge coupling has to be known in the range from μτ F to μT .We can determine the MS gauge coupling g 2 MS = 4πα s self-consistently from the lattice data over a range of scales from about µ = 1/a to about µ = πT by evaluating the squared field strength E = Tr G µν G µν /2 under gradient flow, which yields Here, ⟨E⟩ τF is the expectation value of the squared field strength E at flow time τ F .Then we use the relation between this quantity and the MS coupling, calculated at NLO by Lüscher [25] 20)) as a function of gradient flow scale µF in temperature units.Note that the couplings and scales are defined such that µF = μ.The dotted vertical line is located at µF = μref (Eq.( 21)) and depicts the point from where we move to other scales via the perturbative beta function.Statistical errors are smaller than the linewidth.
Neumann [29]: The MS coupling α s is then obtained by solving the cubic equation.Note that N f is the number of light quark flavors which is zero in the case of this study.This approach becomes unreliable for τ F < a 2 (where a means lattice spacing) due to missing O(a 2 /τ F ) effects.It also has problems when τ F becomes too large, as the matching then occurs in a stronger-coupled regime where the perturbative match of Eq. ( 20) becomes less precise.This makes it difficult to establish α s directly over the full range from μT to μτ F .Therefore the best procedure appears to be to use this direct determination at a single value µ F = 1/ √ 8τ F = μref where both lattice spacing and perturbative effects are under control, and to use the perturbative running of the gauge coupling to obtain it in the desired range.This is the approach we will follow in this work.
Since the matching procedure is conducted based on the continuum theory, it is necessary to extrapolate the flowed coupling constant α flow (µ F ) measured on the lattice (see Table II) to the continuum first.Given that the discretization error of the lattice action is of order a 2 , we adopt an Ansatz linear in a 2 .This extrapolation is shown in the left panel of Fig. 1.It can be seen that the Ansatz describes the lattice data well for all relevant flow times.However, the figure shows that the continuum correction is only small for µ F ≲ 8T , correspondingly restricting the choice for μref .In the right panel of Fig. 1, the continuum-extrapolated (cont.extr.)α flow (µ F ) is shown as a solid green curve.Converting it to α s via Eq.( 20) yields the solid orange curve.We select μref ≡ 2πT (21) and have checked that reducing this value by half does not significantly alter the results.From this point we can move to other scales via the perturbative beta function.
In practice, the evolution is calculated using the function crd.AlphasExact(α s (μ ref ), μref , μτ F , N f = 0, N loop = 5) of the RunDec package [38,39].The resulting coupling is then used to compute Z match via Eq.( 18), which we show in Fig. 2 for different choices of μT and μτ F .The vertical dotted lines enclose all flow times we will use when extrapolating to zero flow (the flow time extrapolation window), meaning that the renormalization procedure changes the bare correlator by at most between −15% and +5%.We can also see that the difference between choices for μτ F is almost negligible, while for μT the difference is simply a constant multiplicative factor.The four curves for Z match shown in this figure, which formally differ from each other at the NNLO level, will be used in the next section to obtain four physical correlators via Eq.( 18).These will be treated on an equal footing in the spectral analysis and therefore contribute to the overall systematic uncertainty.However, as we will see later, the impact of the different choices on the determination of κ B is almost negligible.2. The renormalization factor of the color-magnetic correlator, Z match (Eq.( 18)), as a function of flow scale µF in temperature units (left) or as a function of flow radius √ 8τF in inverse temperature units (right) for the four combinations of μT and μτ F defined in Sec.II.Statistical errors are smaller than the linewidth.The difference between the panels is the inverted x-axis since µF = 1/( √ 8τF).

V. CORRELATOR COMPUTATION
We calculate the color-magnetic correlator G B on the lattice using four different lattice spacings a (cf.Table I) and perform bootstrap resampling to carry out the complete analysis procedure (i.e., continuum extrapolation, renormalization and flow time extrapolation of G B , and in the next section spectral reconstruction) on 1000 bootstrap samples. 6 Note that to suppress discretization errors, we treelevel improve [40] the bare correlator by multiplying it with the ratio of the free correlator obtained in the continuum at leading order G norm B (τ ) [12] (which is same as in the color-electric case) and the one calculated in lattice perturbation theory at leading order at finite flow time [17] This is equivalent to including tree-level lattice-spacing errors in both the correlator and gradient flow proce- 6 For each lattice spacing, 1000 samples for the correlator are created by randomly drawing 10,000 measurements at that lattice spacing with replacement, from which then, for each flow time, the corresponding sample correlator is calculated.The continuum extrapolation yields 1000 continuum correlator samples by drawing one sample correlator from each lattice spacing (without replacement).We represent bootstrap sample distributions using their median ±34th percentile, which is what is shown in all subsequent figures that feature data points with error bars.
dures.To make the subtle features of the data visible, we further normalize it by G norm B (τ ).The quantity of interest is therefore the tree-level-improved ratio For simplicity we drop the trivial superscript "latt" in the following.

A. Continuum extrapolation of the bare correlator
In this section we consider the continuum extrapolation of the bare correlators (Item 2 of Sec.II).The continuum extrapolation of G B (τ, τ F )/G norm B (τ, τ F ) requires an interpolation of the Euclidean temporal separations τ T to the same values across all lattices (we use those present on the finest lattice).This is done using cubic splines with natural boundary condition at the smallest τ T and zero slope at τ T = 0.5 due to the periodic boundary conditions used in the lattice simulation.The continuum extrapolation is carried out using an Ansatz linear in a 2 as we use the standard Wilson action.
In the two panels of Fig. 3 we show the extrapolation at the minimum and maximum flow time [in units of √ 8τ F /τ , cf.Eq. ( 24)] that is later used in the flow-timeto-zero extrapolation.We can see that the linear-in-a 2 fit Ansatz gives a good description for the lattice data at all τ T considered (smaller τ T are not reliably accessible by the gradient flow method, as explained further below).The continuum-extrapolated results are shown in the left panel of Fig. 4. The renormalized correlators obtained via Eq.( 18) are shown in the right panel of Fig. 4 24).The dashed lines and data points at 1/N 2 τ = 0 represent the linear-in-a 2 continuum extrapolation.Statistical errors are smaller than the error bar linewidth.
ameliorate the logarithmic flow time divergence, rendering the remaining flow time dependence of the renormalized color-magnetic correlator linear, which is in close resemblance to the behavior of the color-electric correlator (cf.Refs.[14,17]).

B. Flow-time-to-zero extrapolation of renormalized correlator
Now we perform the τ F → 0 extrapolation of the continuum color-magnetic correlator, which has been renormalized via Eq.( 18) using the matching factors shown in Fig. 2. As detailed in a previous study of the colorelectric field correlators [14,17], there is a small flow time window that is usable for the flow extrapolation, which reads: The lower limit effectively results from the largest lattice spacing used in the continuum extrapolation, as enough flow needs to be applied to sufficiently suppress lattice discretization artifacts.The upper limit results from a perturbative study [41] that allowed the flowed (colorelectric) correlator to deviate from its zero flow counterpart by, at most, 1%.This ensures that the flowed correlator does not stray too far from the physical, zero flow time value.Due to Eq. ( 24), the correlator can only be calculated reliably for Euclidean time separations τ T ≥ 0.25.
In the window given by Eq. ( 24) a linear extrapolation is performed, inspired by the lowest-order small flow time expansion and justified by our data.Such a simple Ansatz turns out to work very well as the right panel of Fig. 4 indicates, where the fits, indicated by the dashed lines, go through all data points in the flow time window.We perform this extrapolation by fitting to the three explicitly shown data points at the minimum, midpoint, and maximum of our flow time window.We do this because there is a large autocorrelation between nearby flow times.Fitting to many nearby flow times is therefore not justified statistically.It is also important that we determine the error in the fit based on the different fit results found in bootstrap resamplings of our data, rather than through naive regression.
In Fig. 5 we show a comparison of the doubleextrapolated color-electric field correlator G E from Ref. [14] and the color-magnetic field correlators G phys.

B
obtained in this study at the same temperature using the same lattices.The figure shows that each correlator is larger at τ T = 0.5 than at small τ T values.This is partly a result of infrared, thermal contributions to the spectral function, and partly because of the intrinsic separation dependence of the vacuum correlation function.
For G E these arise because G E ∝ g 2 which is scale dependent and becomes smaller at short distances-hence the falloff in G E between, say, τ T = 0.35 and τ T = 0.25.For G B this distance dependence arises from g 2 and also from the scale (and hence τ ) dependence of Z match , or equivalently, of c 2 B as we will discuss around Eq. ( 26).This softens the τ dependence of the vacuum The bands depict the statistical ±1σ error around the median value.Right: the same as on the left, but renormalized by multiplying with Z match (Eq.( 18)) for the case μT = μT,NLO, μτ F = μτ F ,NLO (other choices differ only marginally or by a simple rescaling).The dashed lines depict the linear-in-τF/τ 2 flow time extrapolation, which, due to large autocorrelation between subsequent flow times, is based solely on the three explicit data points indicated between 0.25 ≤ √ 8τF/τ ≤ 0.30 (Eq.( 24)).

VI. SPECTRAL RECONSTRUCTION
In this section we invert Eq. ( 5) to reconstruct the spectral function from the renormalized color-magnetic field correlators G phys.  [14]) and color-magnetic (X = B, this work) correlators as a function of Euclidean temporal separation τ in inverse temperature units.For the color-magnetic case four variants are shown that arise from the renormalization due to a priori unknown scales at NLO as detailed in Sec.II.spectral function consisting of two main parts: the ultraviolet (UV) part and the infrared (IR) part.The UV part is known at both leading order (LO) [19] and nextto-leading order (NLO) [19].At LO the MS spectral function matches the color-electric case and reads with C F = (N 2 c − 1)/(2N c ), N c = 3, and g 2 the MS coupling, which should be evaluated at some scale which may be ω dependent.In a previous paper [21] we use this LO spectral function for the UV part, picking a renormalization point which transitions from a fixed value in the IR to μ = 2ω in the UV.However, this approach is not really self-consistent, since it takes into account one source of NLO logarithmic scale dependence, the running of the coupling, without taking into account another source of logarithmic scale dependence, the renormalization of the B-field operators.This operator renormalization also causes a logarithmically scale-dependent shift which is not captured by using an RG-flowed g 2 value in Eq. (25).Therefore, in this work we will only study the case where the UV is modeled by the NLO spectral function.
The expression for the thermal contribution is lengthy, but it turns out to be tiny compared to the vacuum contribution when ω ≫ T . 7At very small frequencies ω < T it is not reliable and in any case small compared to the infrared contribution which the data require us to add.We include the thermal part but have checked numerically that it has only a marginal effect on the spectral reconstruction.
The quantity in round parenthesis in Eq. ( 26) represents the MS spectral function.This should be converted to the spectral function for the physical correlator to account for the renormalization of the magnetic field operators on the Wilson line.This necessitates the factor c 2 B , which plays the same role as Z match played in Sec.II and in particular in Eq. (11).Including this factor, the expression in Eq. ( 26) is formally independent of the renormalization scale μω up to NNLO effects, since the explicit logarithm cancels with the implicit µ dependence from the one-loop running of g 2 (μ ω ) and the µ dependence of c 2 B .Nevertheless, we have to choose some value, and we should attempt to pick a value which will somehow reflect the intrinsic scales present in the spectral function.One such scale is the temperature, and the other is the frequency itself.To estimate the typical thermal scale, we adopt the thermal scale of the dimensionally reduced 3D effective field theory [42], with N c = 3 and N f = 0 yielding μDR ≈ 6.74T .This scale should be used for ω ∼ T .For the relevant scale at large frequency we use μ = ω, and we transition from 7 Unlike the stress tensor spectral function, the color-magnetic and color-electric spectral functions first deviate from their vacuum forms at NLO in the gauge coupling expansion.Structurally, this is because they involve the correlation of a single field strength with another, rather than a product of field strengths, so at the free level the correlation equals the sum of the vacuum-theory correlations summed over periodic copies.The result is that the thermal part first arises at NLO and is therefore small compared to the LO vacuum contribution.
one to the other via To calculate c 2 B we also again need the MS coupling constant, which we reuse from Sec. IV.To take the systematics from higher-loop corrections into account we allow a rescaling of the UV part by the fit parameter K, that is, we choose ϕ UV = Kϕ NLO UV .The IR part manifests a very simple structure based on the infrared asymptotics [43], Now that both IR and UV parts have been established, the missing link is the region between them.Several schemes have been considered in the literature [13-17, 19, 20].In this study we employ three of them that we consider representative.The first one is simply to take the maximum of the two asymptotics [13] ρ max ≡ max(ϕ IR , ϕ UV ). ( The second one makes a smooth switch between the two parts [13] The third one imposes a power-law transition between ϕ IR and ϕ UV [15] The choice for the cuts ω IR and ω UV can be justified based on some physical arguments.Common option of ω IR ranges from T [42] to πT [44].However, according to Ref. [14], ρ max and ρ smax also tend to exhibit a transition around πT .This makes the choice πT redundant so we use ω IR = T .As for ω UV we choose 2πT , which is the standard thermal scale obtained from the perturbation theory in the MS scheme [42].

B
(cf.Fig. 5) as a function of Euclidean time τ for the case μT = μT,NLO, μτ F = μτ F ,NLO (other choices differ only marginally).Different models are slightly offset in τ T for visibility.
With these models we fit the lattice data by minimizing where δG B (τ ) denotes the error of the lattice data and G model (τ ) is calculated using the convolution (Eq.( 5)) replacing ρ B by the above models.The two free parameters of the spectral function models are K and κ B /T 3 .The results of the spectral function model fits are shown in Fig. 6.The left panel highlights the differences between the spectral function model shapes [defined in Eqs.(30) to (32)] in the transition from infrared to ultraviolet asymptotics.The right panel of Fig. 6 shows the residual deviation between the resulting fitted model correlator and the physical color-magnetic correlator data obtained in the previous section.We remark that all the fit procedures in this study are performed on each of the 1000 bootstrap samples separately, with the individual fits using inverse squared statistical errors employed as weights.All models describe the lattice data well according to the residual χ 2 /d.o.f., which is 1.8 ± 0.8 for ρ max , 1.6 ± 0.8 for ρ smax , and 1.2 ± 0.7 for ρ plaw .There is no meaningful difference in the fit quality between the different choices of μT and μτ F .The fit results for the K-factors are K ≈ 0.85 ± 0.01 for all models using μT = μT,LO and K ≈ 1.03 ± 0.01 for all models using μT = μT,NLO .
By aggregating the fit results for κ B /T 3 across different models, we obtain a scatter plot, as shown in Fig. 7.The systematic and statistical errors, combined according to the methodology detailed in [17], yield a confidence interval for κ B /T 3 , which is depicted as a gray band.Including all four combinations for μT and μτ F yields the final estimate In Fig. 8 we compare κ B obtained in this work (black line) and two other calculations, one using the multilevel method [19] and another also using gradient flow with analytical continuation performed at finite flow depth [20].We find that the results are consistent even though they are calculated in different ways.while "TUMQCD '22 (flow * )" refers to Ref. [20], in which gradient flow is the spectral reconstruction is performed at finite flow time, indicated by an asterisk in the legend.
In Fig. 9 we collect all existing determinations for κ E and κ B on the market and plot them logarithmically as a function of g 2 (µ = 2πT ).The quenched data for κ E are taken from Refs.[13][14][15][16]20], while the quenched data for κ B is taken from Refs.[19,20] and this work.The full QCD data for κ E and κ B are taken from Ref. [17] and Ref. [21], respectively.Based on perturbation theory [8], both κ E and κ B are expected to be proportional to g 4 .This is despite the UV behavior of the spectral function scaling as g 2 , as shown in Eq. (25).When we fit all κ B and κ E data using two different Ansätze, one proportional to g 2 and the other proportional to g 4 , we find that the former can describe neither κ E nor κ B and neither quenched nor unquenched data (red curve).But a single fit of form g 4 actually describes all data (quenched and unquenched, κ E and κ B ) within errors (orange curve) with a residual χ 2 /d.o.f.= 0.36.This corroborates the validity of the perturbative prediction, which has also recently been used in an out-of-equilibrium study on heavyquark diffusion [45].

VII. RESULTS AND CONCLUSIONS
In this paper we calculate the mass correction to the heavy-quark momentum-diffusion coefficient using quenched lattice simulations.The correction is extracted from the color-magnetic field correlation functions measured under gradient flow.In our opinion the most interesting part of this paper is Sec.II, in which we elaborate Log-log plot of κE/T 3 and κB/T 3 versus the coupling strength g 2 (μ = 2πT ).The data, from Refs.[13-16, 19, 20] and this work, follow a single trend and turn out to fit well to a function proportional to g 4 , given by κE,B/T 3 = 0.4(0.03)g 4 (orange curve, with the band being the error).The red curve shows an unsuccessful fit of form κE,B/T 3 ∝ g 2 .The fit results are consistent with the perturbative picture that suggests κE,B/T 3 ∝ g 4 [8].
how to take care of the anomalous dimension existing for a magnetic field operator that is further complicated by the matching from the gradient flow scheme to its physical value.By introducing a matching factor Z match , we establish a multistep routine that bridges the gap between the object measured under gradient flow and its physical counterpart.κ B obtained from the physical correlation functions in this work is consistent with previous quenched results calculated in different ways.We also find that κ E,B calculated on the lattice, in both the quenched and the unquenched case, follows a form κ/T 3 = 0.4g 4 (μ = 2πT ) that shows the same scaling with the coupling strength as the perturbative calculation does [8].This remarkably simple expression describes all quenched and unquenched lattice results with a good χ 2 .The methodology developed in this work can be used in the future for a full QCD study at the physical pion mass reaching lower temperatures around the crossover region.
All data from our calculations, presented in the figures of this paper, can be found in Ref. [37].

FIG. 1 .
FIG. 1. Left: gradient flow coupling α flow as a function of squared lattice spacing [or equivalently N −2 τ at fixed temperature T = 1/(aNτ )] at selected flow scales µF in temperature units.The dashed lines correspond to linear fits at fixed flow scale, yielding the continuum-extrapolated data points at N −2 τ = 0. Statistical errors are smaller than the data linewidths.Right: continuum-extrapolated flow coupling α flow and matched MS coupling αs (via Eq. (20)) as a function of gradient flow scale µF in temperature units.Note that the couplings and scales are defined such that µF = μ.The dotted vertical line is located at µF = μref (Eq.(21)) and depicts the point from where we move to other scales via the perturbative beta function.Statistical errors are smaller than the linewidth.

FIG. 5 .
FIG.5.Comparison of renormalized, continuum, zero flow time color-electric (X = E, from Ref.[14]) and color-magnetic (X = B, this work) correlators as a function of Euclidean temporal separation τ in inverse temperature units.For the color-magnetic case four variants are shown that arise from the renormalization due to a priori unknown scales at NLO as detailed in Sec.II.

1 FIG. 7 .
FIG. 7. Coefficient κB of the finite-mass correction to heavyquark momentum diffusion from each model.The gray band shows the confidence interval that ranges from the 34th to 68th percentile of the total distribution of the fit results of all bootstrap samples for all models (see Sec. V).

FIG. 8 .
FIG.8.Comparison of quenched κB/T 3 obtained in this work (including the statistical and systematic uncertainty from all four choices of μT and μτ F , as well as the various spectral function models of Sec.VI) and recent literature."Banerjee '22 (ML)" refers to the multilevel method work of Ref.[19], while "TUMQCD '22 (flow * )" refers to Ref.[20], in which gradient flow is the spectral reconstruction is performed at finite flow time, indicated by an asterisk in the legend.
s via Eq.20 at µ F = μref , then run via 5-loop β-fct.α s via Eq.20 α flow (cont.extr.) α G B correlator in comparison to G E .FIG.4.Left: bare continuum-extrapolated color-magnetic correlator GB, tree-level-improved and normalized to its free counterpart G norm , as a function of normalized flow time τF/τ for Euclidean temporal separations τ in units of inverse temperature.
Left: fitted model spectral functions ρB function of frequency ω in temperature units for the case μT = μT,NLO, μτ F = μτ F ,NLO (other choices differ only marginally).Statistical errors are not shown to reduce clutter; the intent of the figure is to show the differences in spectral function model shapes.The y-axis is scaled to yield κB/T 3 for ω → 0. Right: ratio of fitted model correlators G model to renormalized, continuum, zero flow time correlator data G phys.