Vacuum correlators at short distances from lattice QCD

Non-perturbatively computing the hadronic vacuum polarization at large photon virtualities and making contact with perturbation theory enables a precision determination of the electromagnetic coupling at the $Z$ pole, which enters global electroweak fits. In order to achieve this goal ab initio using lattice QCD, one faces the challenge that, at the short distances which dominate the observable, discretization errors are hard to control. Here we address challenges of this type with the help of static screening correlators in the high-temperature phase of QCD, yet without incurring any bias. The idea is motivated by the observations that (a) the cost of high-temperature simulations is typically much lower than their vacuum counterpart, and (b) at distances $x_3$ far below the inverse temperature $1/T$, the operator-product expansion guarantees the thermal correlator of two local currents to deviate from the vacuum correlator by a relative amount that is power-suppressed in $(x_3\:T)$. The method is first investigated in lattice perturbation theory, where we point out the appearance of an O$(a^2 \log(1/a))$ lattice artifact in the vacuum polarization with a prefactor that we calculate. It is then applied to non-perturbative lattice QCD data with two dynamical flavors of quarks. Our lattice spacings range down to 0.049 fm for the vacuum simulations and down to 0.033 fm for the simulations performed at a temperature of 250 MeV.

Non-perturbatively computing the hadronic vacuum polarization at large photon virtualities and making contact with perturbation theory enables a precision determination of the electromagnetic coupling at the Z pole, which enters global electroweak fits. In order to achieve this goal ab initio using lattice QCD, one faces the challenge that, at the short distances which dominate the observable, discretization errors are hard to control. Here we address challenges of this type with the help of static screening correlators in the high-temperature phase of QCD, yet without incurring any bias. The idea is motivated by the observations that (a) the cost of high-temperature simulations is typically much lower than their vacuum counterpart, and (b) at distances x 3 far below the inverse temperature 1/T , the operator-product expansion guarantees the thermal correlator of two local currents to deviate from the vacuum correlator by a relative amount that is power-suppressed in (x 3 T ). The method is first investigated in lattice perturbation theory, where we point out the appearance of an O(a 2 log(1/a)) lattice artifact in the vacuum polarization with a prefactor that we calculate. It is then applied to non-perturbative lattice QCD data with two dynamical flavors of quarks. Our lattice spacings range down to 0.049 fm for the vacuum simulations and down to 0.033 fm for the simulations performed at a temperature of 250 MeV.

I. INTRODUCTION
Many phenomenologically interesting observables are defined in terms of QCD vacuum correlators involving two or more local fields integrated over their Euclidean positions. For example, the hadronic vacuum polarization function Π(Q 2 ), which determines the leading hadronic contribution to the running of the electromagnetic coupling and the muon anomalous magnetic moment (g − 2) µ , the π 0 → γ * γ * transition form factor of the pion and the hadronic light-by-light contribution to (g − 2) µ are expressed in such a way. Thus in many cases, vacuum correlators represent crucial input for precision tests of the Standard Model. Lattice QCD provides ab initio determinations of these vacuum correlators; see e.g. Refs. [1][2][3][4] for the applications above.
However, these integrated quantities contain contributions corresponding to the fields being close together, a regime which can lead to large cutoff effects. The standard tool to investigate the asymptotic approach to the continuum limit of correlation functions is Symanzik's effective field theory [5][6][7]. A complication for the aforementioned observables is that the on-shell improvement programme is not sufficient to guarantee rapid convergence toward the continuum limit. In this work, we propose to compute the short-distance contribution to vacuum correlators by making use of static screening correlators from QCD at finite temperature, which can significantly reduce the cost of obtaining a robust continuum limit. As the short-distance contribution has little sensitivity to the temperature, the bulk of this contribution can be computed using particularly small lattice spacings in the hightemperature phase of QCD, where the cost of the simulations is much reduced, and only a small remainder needs to be computed using the vacuum ensembles. Just as importantly, the cutoff effects on the remainder can be arranged to be parametrically smaller than those of the observable computed on thermal gauge ensembles. Furthermore, we note that in certain cases, there is a logarithmic enhancement of the cutoff effects on the short-distance contribution at leading order in perturbation theory, in contrast to the modification of cutoff effects by logarithms affecting on-shell correlators, which only appears beyond the free-field theory level [8]. The longer-distance contribution involves the currents at physical separations, at which on-shell improvement can safely be applied directly to the vacuum correlators.
The fact that the leading thermal effect on the correlator at short distances x 3 is suppressed by several powers of (x 3 T ) allows for a strategy to compute correlators at extremely high momentum scales |Q|. We make a fairly concrete proposal in this direction in Section V for how to compute the hadronic contribution to the running of the electroweak coupling constants up to the Z-boson mass. The basic idea is to compute this contribution by increasing the momentum scale by factors of two, always using thermal QCD ensembles with |Q|/T sufficiently large that the thermal effects are small corrections computable to a systematically improvable accuracy. The idea thus has common aspects with the 'step-scaling' idea introduced in the lattice field theory context in [9], and also with the application in heavy-quark physics presented in Ref. [10].
In the following section, we outline the strategy for computing short-distance observables using auxiliary finite-temperature ensembles, and provide parametric estimates for the optimal choice of lattice parameters. We examine in Section III the case of the vector current correlator in the free theory, which suggests that the thermal effects are guaranteed to be small at sufficiently small separations of the currents. This is confirmed by the operator product expansion, carried out at leading and next-to-leading order in Appendix A. In Section IV, we test the strategy on vector current correlators with N f = 2 Wilson fermions with thermal ensembles with a temperature T = 250 MeV, where we reach lattice spacings of a ≈ 0.033fm. This provides a more controlled continuum limit of the short-distance contribution to the vacuum polarization at a much reduced cost. Finally, Section V summarizes our findings and describes the idea to compute the running of electroweak couplings to very high energy using sequences of ensembles of growing temperature. The concrete setup suggested is tested in the free-theory context in Appendix E. The other appendices contain technical details of the analytic calculations.

II. DEFINITIONS & GENERAL IDEA
To be specific, in this study we concentrate on two observables which are defined in terms of the integral over a Euclidean correlator weighted by a known kernel and are closely related to the hadronic vacuum polarization function Π(Q 2 ). The Adler function [11] parametrizes the running of the hadronic contribution to the electromagnetic coupling at spacelike q 2 = −Q 2 < 0, and its derivative at Q 2 = 0 determines the anomalous magnetic moment of a lepton in the limit of vanishing lepton mass, m l . Both of these quantities receive contributions from all non-zero time-separations of the current correlator where the (continuum) electromagnetic current is defined as and Q f is the electric charge of quark flavour f = u, d, s, . . . and the matrices γ µ satisfy the Euclidean Dirac algebra {γ µ , γ ν } = 2δ µν . The kernels for both the Adler function and the anomalous lepton magnetic moment (g − 2) µ in the time-momentum representation [11] coincide with the fourth moment for small enough x 0 . In fact, in the case of (g − 2) µ , the kernel agrees with the fourth moment at the percent level up to distances of about 0.5 fm. Therefore, both the qualitative and quantitative results for the fourth moment are relevant for the controlled determination of the short-distance contribution to the hadronic vacuum polarization in the muon anomalous magnetic moment. One may wonder whether lattice QCD estimates of these physical quantities suffer from uncontrolled systematic effects arising from small separations between the currents, even if this contribution itself is suppressed by the short-distance behaviour of the kernel. Indeed, for a lattice estimator of the current correlator which has the expansion in the lattice spacing a given by power counting suggests that its cutoff effects become parametrically large as x 0 becomes small [12] a n G n (x 0 , a) = const. × (a/x 0 ) n G(x 0 ) + . . . , up to logarithmic corrections [8], and where we assume the continuum limit G(x 0 ) exists after proper renormalization, if required. Even if the bulk of the lattice artifacts does not necessarily arise from the short-distance contribution, the breakdown of the Symanzik expansion inevitably leads to scaling violations in the continuum limit which can be of practical concern, especially given the subpercent precision aimed at in the context of (g−2) µ . integrand, a = 0.05 fm  This situation is similar to the typical window-problem encountered in lattice QCD where the appearance of an external scale, such as Q 2 , needs to be accommodated within the ultraviolet and infra-red cutoffs imposed by the lattice spacing a and lattice size L, In renormalization problems, one has the freedom to remove one of these restrictions by linking the external scale to the physical volume, which eliminates one constraint of the window and allows simulations to proceed with tractable problem sizes. For hadronic observables, where the physical volume must remain large, we may however choose to compute an observable, or part of it, in a simulation with different physical parameters provided that we properly account for the correction. In particular, our strategy proposes to use the static screening correlator at finite temperature T , which is a function of the spatial separation x 3 of the currents and depends on the temperature T . We choose the latter to be on the order of the QCD scale, in the chirally-restored phase. We define the contribution up to t to the integral appearing in Eq. (3) for the vacuum and thermal correlators, together with lattice estimators I(t, a) and its thermal counterpart I th (t, a) defined precisely in the following subsection. Our strategy is based on the idea that the quantities I(t) and I th (t) are in some sense very similar. The operator-product expansion (OPE), which is presented in detail in Appendix A, can be invoked to make this statement precise for tT 1: the difference of the two quantities is suppressed by (tT ) 3 relative to the quantities themselves. It is instructive to compare the thermal and the vacuum correlators in a representative lattice QCD calculation. The left panel of Figure 1 depicts the integrand of Eq. (10) with t = 0.2 fm for the vacuum (open) and thermal (filled) squares at fixed lattice spacing a ≈ 0.05 fm, which illustrates that the thermal effects are indeed suppressed for these distances in the N f = 2 theory when T = 250 MeV, corresponding to tT = 0.25. The analogous integrands for the Adler function at the large virtuality of Q = 2.36 GeV are shown as well, which illustrate how it is also dominated by the correlation function at short distances. The right-hand panel shows the corresponding integrals as a function of the lattice spacing, which illustrates that more than 95% of the signal is accounted for by the thermal observable. The benefit of using I th (t) as a proxy for I(t) is that, in the case illustrated in Figure 1, the finite-temperature ensemble has a factor eight fewer lattice sites than its vacuum counterpart due to its shorter time extent, which naively allows a span of a factor of 8 1/4 ≈ 1.68 in the lattice spacing to be achieved for the thermal observable before the cost of obtaining the latter becomes comparable to the vacuum calculation.
Thus for a suitable choice of t 1/T , we expect the bulk of the short-distance contribution to be given by the thermal component I th (t), whose continuum limit can be obtained accurately thanks to the smaller lattice spacings accessible at finite temperature. This suggests an improved estimator for the vacuum observable where the first term on the right-hand side is the continuum estimate of the thermal observable, obtained using particularly fine lattice spacings available at high temperature. The remainder in brackets is small and of the form const. × t 5 (1 + O(t)), where the constant is dominated by a momentum scale on the order of temperature. It is worth recording the parametric size of the cutoff effects on both terms: for the first term, the cutoff effects are of order a 2 , while for the remainder, they are of order a 2 (T t) 3 in the O(a)-improved theory 1 . Thus, as long as the ratio of the lattice spacings used in the thermal theory to those used in the vacuum does not become as small as (tT ) 3/2 , which is the regime we have in mind, the cutoff effects on I th (t, a) are parametrically larger than the effects on the remainder. The upshot is that if I(t) is obtained as the continuum limit of I(t, a), the dominant part of the systematic error associated with cutoff effects comes from obtaining I th (t), where one profits from being able to reach lattice spacings below 0.05 fm in the chirally-restored phase at a moderate computational cost. A further aspect which is specific to Wilson fermions is that the on-shell improvement of the vector current via the derivative of the tensor current [13] only contributes a term of order am q in the chiral restored phase of QCD; such terms are of a size comparable to the O(a 2 ) terms for the lattice spacings employed in Section IV. This represents a further advantage of obtaining the bulk of I(t, a) from the chiral restored phase.
Since the main focus is then on obtaining I th (t), it is worth studying the approach to the continuum for this short-distance quantity in lattice perturbation theory. The leading-order calculation is presented in Section III. It turns out that a logarithmic enhancement of the O(a 2 ) cutoff effects arises, with a calculable coefficient which applies both to I th (t) and I(t). This enhancement appears with a positive (unit) power of the logarithm, unlike the known logarithmic dependence on the lattice spacing due to the running coupling, which appears first at one-loop level [8]. By contrast, the O(a 2 ) cutoff effect enhanced by the factor log(1/a) cancels out in the improved observable (11).
In order to fully control the short-distance thermal contribution, or to reach very high momenta in the hadronic vacuum polarization, it may be necessary to iterate the procedure of Eq. (11) using a series of higher temperatures to compute short-distance contributions. We return to this question in Section V. Although other options are certainly available, it is particularly convenient to use the temperature as a control parameter to compute the short-distance contribution, since we can use existing knowledge about high-temperature correlators and apply well-understood theoretical tools like the operator product expansion.

A. Definitions of lattice observables
In order to set up the notation for the following sections, we define here the lattice observables for the theory of N f = 2 Wilson fermions. In this work we investigate the isovector vector current correlator, which consists of a single Wick contraction. This correlator makes the dominant contribution to the hadronic vacuum polarization in the muon g − 2. In the vacuum case, we formulate the correlator at vanishing spatial momentum as a function of Euclidean time, where the bare local vector current is defined as with Ψ = (u, d), and the exactly-conserved vector current is In contrast, the static screening correlator at finite temperature is measured along a spatial direction, We investigate two observables which are related to the Adler function Eq. (1), and the short-distance part of the fourth moment of the correlator Eq. (10), which is proportional to the derivative of the vacuum polarization at zero virtuality. The short-distance contribution up to t of the fourth moment is We have used an integration rule consistent with the improvement of the theory, e.g. the trapezoidal rule. For large Q 2 , the Adler function is dominated by the short-distance contribution to the integral, and we define a lattice observable which is the integral up to half the spatial extent L/2, where again we have implemented the trapezoidal rule, and K(x 0 , Q 2 ) is defined in Eq. (2). The thermal observables are defined analogously using the static screening correlator of Eq. (15). The improved estimators are then defined via Eq. (11) where the first term on the right-hand side is obtained by taking the continuum limit using thermal ensembles with the available finer lattice spacings.

B. The Symanzik expansion and enhanced lattice artifacts
As discussed at the beginning of the section (see Eqs. (6-7)), severe lattice artifacts appear in the correlation function at short distances. Here, we demonstrate that integrating over the correlation function with a kernel suppressing the short distances sufficiently so as to yield a finite continuum limit can result in a parametric enhancement of lattice artifacts even at leading order in the perturbative expansion.
The Symanzik continuum effective theory can be used to represent the correlation function on the lattice for x 0 > 0 by considering all irrelevant counterterms of the action and local operators with the correct dimension and consistent with the symmetries of the lattice theory. For example, the lattice artifacts of Eq. (7) can be expressed as a sum over the matrix elements containing the counterterms [8] with coefficients which depend on the (scheme-independent) one-loop anomalous dimension of the counterterm γ i 0 (b 0 is the universal one-loop coefficient of the QCD beta function) andc i is the matching coefficient between the Symanzik continuum effective theory and the lattice theory. The matrix element C i n is renormalization-group invariant as the scaledependence of the counterterm has been factored out, which gives rise to a logarithmic dependence on the lattice spacing through the running coupling.
In addition, however, the integral of the correlator from short-distances results in a logarithmically-enhanced lattice artifact. In Eq. (7), the matrix element of any one of the leading O(a 2 ) counterterms must have by power counting the short-distance singularity This is more singular than the leading continuum correlator, and gives rise to a logarithmic enhancement of the O(a 2 ) lattice artifacts in the x 4 0 moment, in particular when inserted in the summation of Eq. (16), using the harmonic number formula. Thus, assuming the correlator is O(a)-improved, the leading O(a 2 ) lattice artifacts of the integrated quantity are parametrically enhanced due to the logarithm appearing with the (positive) unit power. It is also worth noting that all higher terms in the Symanzik expansion of G(x 0 , a) contribute at O(a 2 ) after summing over short distances. In particular, even if one improved the onshell correlator so that it contained no O(a 2 ) artifacts, the quantity I(x 0 , a) would contain remnant O(a 2 ) lattice artifacts.
The full form of the lattice artifacts given by Eq. (18) suggests that the coefficient of the logarithmic term could be determined by computing all of the matching, or improvement, coefficients while the most singular behaviour of the coefficient function is computable in continuum perturbation theory, when x 0 Λ 1. In the following section, the coefficient of the logarithmically-enhanced term is computed at leading order in lattice perturbation theory.

III. ANALYSIS IN LEADING ORDER OF LATTICE PERTURBATION THEORY
As a first test of the idea to use thermal gauge ensembles to better control the shortdistance behaviour of QCD correlators, we apply it in the framework of leading-order lattice perturbation theory, i.e. in the theory of non-interacting quarks. In the limit of short distances, the QCD correlation functions are well approximated by their perturbative values, and we therefore expect the free theory to provide valuable insight on the general applicability of the method. Moreover, the continuum values being known in this context, quantitative statements can be made about the accuracy of the continuum limit obtained with the improved estimators in Eq. (11) as compared to extrapolating directly the vacuum observables.
One delicate point in the study of the integrated observables defined in Section II A is the presence, already at the level of the free theory, of cutoff effects which depend logarithmically on the lattice spacing. In the method we propose, the understanding of these effects is important in view of obtaining an accurate continuum extrapolation of the thermal observables. On the other hand, the improved vacuum observables defined as in Eq. (11) are free of the logarithm predicted by the leading-order calculation, which cancel out in the subtraction between thermal and vacuum quantities.

A. The vector correlators in the massless theory: lattice formulation
We consider the theory of non-interacting massless Wilson quarks, defined on a lattice with infinite spatial volume. Following Section II A, we denote the lattice vacuum and thermal correlation functions by G µν (x 0 ) and G th µν (x 3 ) (Eqs. (12) and (15), recalling that Z V = 1 in the non-interacting case). Explicit expressions of the free correlators are given in Appendix D for the theory with N c colors. Here we fix N c = 3, as appropriate for QCD. We concentrate on the observables I(t) and D(Q 2 ), as defined in Eqs. (16) and (17).
The analysis is performed at a set of realistic lattice spacings which correspond to those available in our non-perturbative study of N f = 2 QCD (see Table IV). In the free massless Lagrangian there are no bare parameters to be tuned in order to approach the continuum on a line of constant physics, given that the mass parameter is only multiplicatively renormalized in this case. The fact that the non-interacting massless theory is scale invariant gives us the freedom to assign to the temperature a value of our choice. In the thermal case, the lattice spacing and the temperature are related by T = 1/(aN t ), where N t = L 0 /a is the number of lattice points in the Euclidean-time direction. As in the case of the interacting ensembles, we consider four finite-temperature lattices, with N t = 12, 16, 20, 24, and we Free massless quarks vacuum finite temperature assign to each of them the physical temperature T = 246.25 MeV, which corresponds to fixing the lattice extent in the compact direction to aN t = 1/T = 0.8 fm. This results in the set of lattice spacings a ≈ {0.07, 0.05, 0.04, 0.03} fm. Vacuum equivalents of these thermal systems are obtained by assigning the corresponding physical value of the lattice spacing to a lattice with N t = ∞ (zero temperature). For simplicity, with an abuse of notation we will sometimes identify the vacuum lattices by the value of N t of their thermal counterpart (as for example in Figure 4). We analyze I(t) for t = 1/(4T ) = 0.2 fm and t = 1/(2T ) = 0.4 fm, and for the Adler function we consider the two virtualities Q = 3πT = 2.32 GeV and Q = πT /2 = 387 MeV. We are mostly interested in the more short-distance-dominated cases t = 0.2 fm and Q = 2.32 GeV, for which the method proposed in this paper proves to be very effective. Similar values of t and Q are used in the analysis of the N f = 2 QCD data presented in Section IV. The more infrared scales t = 0.4 fm and Q = 387 MeV are only considered in the free-theory analysis, and they mostly serve as a comparison point. Figure 2 shows the integrand of the fourth-moment observable x 4 µ G(x µ ) for all thermal lattices and their vacuum counterparts, up to distances of 0.4 fm. As expected, up to around 0.2 fm the difference between the vacuum and thermal cases is hardly noticeable. Also, at these short distances the cutoff effects are more important than at larger separations, as can be observed by comparing with the continuum values, also displayed in the plot. These two features motivate the use of the improved observables defined as in Eq. (11) in order to achieve a better control at short distances with the aid of fine thermal lattices.
Before undertaking the analysis of the observables I(t) and D(Q 2 ) with the method proposed in this paper, we investigate the emergence of logarithmic cutoff effects and compute their form explicitly.

B. A short-distance O(a 2 log(1/a)) cutoff effect
Already at the free-theory level, cutoff effects of the formc a 2 log(1/a) are present in the observables I(t) and D(Q 2 ), and analogously in the corresponding thermal quantities I th (t) and D th (Q 2 ). In the following we compute the coefficientsc I andc D for the specific discretization of the correlation functions used in this work.
To begin with, we focus on I(t) and we analyze it in the limit a → 0. In the vacuum and in infinite volume there is no difference between projecting to zero momentum in the directions (x 1 , x 2 , x 3 ), as in Eq. (12), or in the directions (x 0 , x 1 , x 2 ). Here, as in Appendix D, we choose the second option, as it makes the analogy with the thermal screening correlator very clear. As a consequence, we will use the vector notation p ≡ (p 0 , p 1 , p 2 ). The thermal version of the following equations is obtained by replacing the integral over p 0 with a sum over fermionic Matsubara modes, as described in Appendix D, and the coefficientsc I and c D are the same in the vacuum and in the thermal case. In the limit a → 0, the observable I(t) can be expanded as follows 2 , where p ≡ |p| andf n,m (p) are dimensionless functions of the orientation of the vector p (p ≡ p/p). The expressions off 00 ,f 2,0 andf 2,1 can be found in Appendix D. A generic term of the expansion within square brackets in (22) can be expressed as a n |x 3 | m p n+mf n,m (p) , with n = 2, 4, . . . and m ≥ 0. The integration over Upon integration over the Brillouin zone, the first term on the right-hand-side of Eq. (24) (not exponentially suppressed in p) can introduce a logarithmic dependence on a. In fact, based on dimensional analysis, we observe that the term is proportional to a 2 log(1/a) for n = 2, and to a 2 for any n > 2. As a consequence, the O(a 2 ) contribution to I(t) cannot be computed exactly by truncating the expansion in square brackets in Eq. (22) at a finite order. Having identified the sources of logarithmic contributions, we can compute the coefficientc I by using the log-derivativẽ The derivative of the triple integral can be computed by applying the following equation Using the expressions off 2,0 andf 2,1 given in Appendix D, we find Moving now to the Adler function D(Q 2 ), we observe that in the short-distance limit its integrand is proportional to that of I(t) As we saw in the above computation, the logarithmic cutoff effect comes from the contribution around x 3 = 0 to the integral over x 3 (see Eq. (24)), from which we conclude thatc The values ofc I andc D given in Eqs. (28) and (30) can be compared to what is obtained via fits to the lattice data. With this goal in mind, we consider a set of lattices with N t = 16, 20, 24, 48, 60, 120, whose temperature is fixed to T = 246.25 MeV. These include "realistic" lattices with N t = 16, 20, 24, whose lattice spacings are very close to those of the ensembles presented in Section IV, and three extremely fine lattices with N t = 48, 60, 120. We include the latter in order to have a better control on the continuum extrapolation and more flexibility with respect to the number of fit parameters. We consider the functional forms (a) c 0 + a 2 [c 2 +c log(1/(T a))] and observe that the case (c) leads to the best agreement with the expected value ofc. In all cases, the agreement between c 0 and the known continuum value is very good. All results are reported in Table I.

C. Continuum limit of the thermal observables
As a first step toward improved vacuum observables defined as in Eq. (11), we compute continuum estimates of the thermal quantities I th (t) and D th (Q 2 ). We consider a set of Table I: Parameters of the continuum extrapolation of I th (t) and D th (Q 2 ) with the functional forms (a), (b), (c), (d) of Eq. (31). The known continuum values I th (t) and D th (Q 2 ) are also reported, together with the prefactors of the logarithmic termc I (Eq. (28)) andc D (Eq. (30)). For this analysis an extended set of thermal lattices is used, with N t = 16, 20, 24, 48, 60, 120 and whose temperature is fixed to T = 246. 25 MeV. For all fit forms, c 0 is in good agreement with the corresponding continuum value. In the more infrared cases t = 0.4 fm and Q = 387 MeV the ansatz (a) provides a good estimate ofc, and the value of this coefficient is quite stable with respect to introducing higher powers of a in the fit ansatz. Instead for t = 0.2 fm and Q = 2.32 GeV the role of higher-order discretization effects is important to obtain a relatively accurate estimate ofc. In particular, among the fit forms analyzed here, the ansatz (c) provides the most accurate value.
Free massless quarks  four lattices with N t = 12, 16, 20, 24, whose temperature is set to 246.25 MeV and whose lattice spacings are very similar to the ones of the N f = 2 thermal ensembles, as discussed in Section III A. In all practical cases, we find it expedient to exclude the coarser lattice with N t = 12 from the continuum extrapolation. We observe that a careful treatment of the logarithmic cutoff effects can significantly improve the accuracy of the continuum limit. In order to illustrate this point, we compare the outcome of two different approaches. The first is to simply ignore the presence of logarithmic cutoff effects and to fit the lattice data with polynomials in a 2 . The second is to subtract the logarithmic contributions making use of the known prefactorsc I andc D (see Section III B) before fitting polynomially in a 2 . This second approach proves to be very effective, however it is somewhat specific to the free theory. In the interacting case, it remains to be seen whether the a 2 log(1/a) term receives significant, non-analytic in a corrections. Another possible strategy is to fit the lattice data with the ansatz c 0 +c 2 a 2 +ca 2 log(1/(T a)).   Table II.
In this case we find that the accuracy of the continuum estimates is comparable with the results obtained by subtracting the logarithmic term and then fitting linearly in a 2 , which are given in Table II. In the left panel of Figure 3 two sets of lattice data are shown, one corresponds to the plain thermal observable I th (t) (open symbols) and the other represents I th (t) −c I a 2 log(1/(T a)) (filled symbols). The value of the coefficientc I is given in Eq. (28). Four different continuum estimates are obtained by fitting these data sets with two polynomial forms, one linear in a 2 and one quadratic in the same variable. The discrepancy between these estimates and the known continuum value is reported on the left-hand side of Table II. The accuracy of the continuum limit is significantly improved by the subtraction of the logarithmic term and the inclusion of the O(a 4 ) term also plays an important role. For example, for t = 0.2 fm a naive polynomial fit of the lattice data gives a discrepancy with the correct continuum value of around 2%, which is reduced to 0.2% by subtracting the logarithmic cutoff effects and fitting the resulting lattice points with a second-degree polynomial in a 2 . As final continuum Table II: Accuracy of the continuum limit of I th (t) (left) and D th (Q 2 ) (right) expressed in terms of the relative difference between the continuum estimate c 0 and the known continuum value. Three thermal lattices are used for the extrapolation, with N t = 16, 20, 24 and whose temperature is set to T = 246.25 MeV. The label 'plain' refers to a continuum estimate obtained by fitting the plain lattice observables I th (t) and D th (Q 2 ), while in the case 'subtr.' the logarithmic cutoff effects are subtracted prior to performing the fit as follows, I th (t) −c I a 2 log(1/(T a)), D th (Q 2 ) − c D a 2 log (1/(T a)). The analytic values of the coefficientsc I andc D are given in Eqs. (28) and (30) respectively. The lattice data and the fit curves are shown in Figure 3.
estimates we choose the most accurate results A similar analysis of the thermal Adler function D th (Q 2 ) can be found on the right panel of Figure 3 and on the right-hand side of Table II. Also for this observable the gain in accuracy due to subtracting the logarithmic cutoff effects is considerable. For example, for Q = 2.32 GeV the continuum estimate obtained by fitting D th (Q 2 ) with a seconddegree polynomial in a 2 differs from the correct continuum value by 0.7%, while fitting D th (Q 2 ) −c D a 2 log(1/(T a)) with the same ansatz reduces the discrepancy to 0.06%. The value ofc D is given in Eq. (30). As final continuum estimates of the thermal Adler function we choose the most accurate values     Table III.

D. Continuum limit of the improved vacuum observables
With the continuum estimates of Eqs. (32) and (33), we build the improved vacuum observablesÎ As in the case of the N f = 2 QCD ensembles, we evaluate the observables on two zerotemperature lattices, whose lattice spacings a ≈ {0.07, 0.05} fm are the same as those of the N t = 12, 16 thermal lattices. We obtain continuum estimates by fitting linearly in a 2 the lattice observables I(t), D(Q 2 ) and their improved versions defined in Eqs. (34) and (35). The lattice data and the fit curves are shown in Figure 4, while the accuracy of the resulting continuum estimates, given in terms of the relative difference with the correct continuum value, is reported in Table III. For the case of I(t = 0.2 fm), one order of magnitude in accuracy is gained by using the improved lattice observables introduced in this paper. In all cases considered here, the advantage of using the improved observablesÎ(t),D(Q 2 ) is rather clear as far as the accuracy of the resulting continuum estimates is concerned.
For the more ultraviolet scales t = 0.2 fm and Q = 2.32 GeV there is the further benefit of a significant reduction of the cutoff effects at finite lattice spacing, as compared to the plain lattice observables I(t), D(Q 2 ). As a final remark, we repeat that the logarithmic discretization effects proportional to a 2 log(1/a) cancel inÎ(t) andD(Q 2 ), due to the subtraction between the vacuum and thermal lattice observables.

IV. NON-PERTURBATIVE TEST IN N f = 2 QCD
In this section, we perform a non-perturbative numerical study of the same observables investigated in the free theory in the previous section, namely the (truncated) fourth moment of the current correlator and the Adler function at large virtuality. We make use of the vacuum CLS ensembles with N f = 2 non-perturbatively O(a)-improved Wilson fermions and the Wilson gauge action with two lattice spacings of a ≈ 0.049 fm and a = 0.0658 fm. Our study is performed at a fixed, common mass of the up and down quarks. For the (zero-temperature) pion mass, we quote the values 268(3) MeV and 269 (3) MeV respectively for ensembles F7 and O7 [14]. The improved estimators were computed using ensembles on the same line of constant physics set by the physical volume L and quark mass m q with a temperature of T = 250 MeV.
The aspect ratio for the finite-temperature ensembles was set to L 0 /L = 1/4, and for the vacuum ensembles to L 0 /L = 2. The thermal ensembles were recently used in a study of the photon emissivity of the quark-gluon plasma [15]. Two of them have common bare parameters with the vacuum ensembles, while two additional ensembles with lattice spacings down to a 0.033 fm allow the continuum limit of the thermal observable to be obtained with reduced uncertainty. Further details on the ensembles are collected in Table IV. 3 The scale was set for the F7 ensemble taking the lattice spacing from ref. [16], and assuming a perfect line of constant physics with a ratio of lattice spacings of 3/4 between O7 and F7. The assumed ratio of lattice spacings is consistent at the one-sigma level with the values of the lattice spacings given in [16], as well as with those of Ref. [14], in which they are quoted with a 0.9% precision. Note that in contrast to the free massless theory, in the present case the current is not fully O(a)-improved as the improvement coefficients are not known non-perturbatively. Nevertheless, at high temperature, the O(a) discretization effects due to the missing current improvement terms should be proportional to the quark mass O(am q ), owing to the restoration of chiral symmetry in the massless theory [17,18].
The integrands for the two observables considered here are displayed in Figure 1 for the vacuum and thermal O7 ensembles. For both observables O = I, D, we employ linear or quadratic fit ansätze and, for the thermal observable, additionally an ansatz where the logarithm is included using the coefficient determined at leading order in the previous section  [19] using the Schrödinger functional coupling computed in Ref. [20].   For the thermal observables we also use the leading-order result of the previous section to implement an additive perturbative improvement according to We now discuss these two observables in turn.  I(t, a), while the filled ones represent the estimator I(t, a) of Eq. (11). Both data sets are fitted linearly either in a or in a 2 . We also estimated the observable I(t) using the perturbative five-loop vacuum spectral function, following Refs. [21,22], depicted with the red point.
A. The short-distance contribution to Π (Q 2 = 0) First we examine the estimate of the continuum limit of the truncated fourth moment of the thermal correlator Eq. (10), which is required to compute the improved estimator Eq. (11). In the left panel of Figure 5, the integral up to t = L 0 /4 = 0.1974 fm is shown as a function of the lattice spacing, for the thermal observable. While the (Q+log) fit provides a satisfactory description of the data when the coarsest lattice spacing is omitted, the perturbatively-improved observable has a much flatter behaviour toward the continuum. The fit results are given in Table IV. The continuum estimate from the leading-order improved extrapolation is used to define the improved estimator of Eq. (11), which is shown in the right panel of Figure 5. The original data are displayed as open symbols; they exhibit a large cutoff effect. For illustration, two fit ansätze are employed, purely linear or purely quadratic in the lattice spacing. The latter may seem more plausible here, given the short-distance nature of the observable. Nevertheless, the severity of the cutoff effect leads to an unsatisfactory control of the continuum limit using only vacuum correlators with the available state-of-the-art lattice spacings. On the other hand, the thermal-improved estimator depicted with filled symbols shows an almost flat continuum limit, which suggests the subtraction of the thermal contribution also removes a significant amount of the cutoff effects, as expected. In this case, the continuum result is much less sensitive to the choice of fit ansatz for taking the continuum limit.
In order to quote a continuum estimate for the thermal-improved observable for illustration, we choose to use the continuum estimate of the thermal observable obtained with the LO-improvement, and for the correction the mean of the continuum estimates obtained with the linear and quadratic fits to arrive at The second, systematic error associated with taking the continuum limit is estimated as the quadrature sum of (a) the difference between the (Q+log) extrapolation of I th and the (Q) extrapolation ofI th , and (b) half the difference between the linear and quadratic continuum fits of the correction term. Finally, we compute the same observable in N f = 2 massless perturbation theory based on the spectral representation [11] where ρ(ω 2 ) is the five-loop massless vacuum spectral function 4 [21,22]. We use the renormalization scale µ = 2.4 GeV and take Λ MS from the FLAG report [23]. The result obtained is I pert (t = 0.1974fm) = 1.059(−6)(1) × 10 −3 fm 2 , which is depicted with the filled red point to the left in Figure 5. The errors are estimated using the uncertainty in Λ (2) MS . We find reasonable agreement between our final estimate and perturbation theory, even though we have not investigated the systematic effects associated with finite quark masses and residual non-perturbative effects, which might become relevant at the quoted level of precision. Also, we have not included the scale-setting uncertainty in Eq. (40); the relative scale-setting uncertainty of I(t) is mainly that of t 2 /fm 2 , i.e. about 2%, whereas I(t)/t 2 only depends weakly on t around t = 0.2 fm.

B. The Adler function at large Q 2
The continuum limit of the thermal Adler function for a large value of the virtuality Q = 3πT = 2.36 GeV is shown in Figure 6. As in the previous case, the leading-order improved observable shows better scaling to the continuum limit, when the coarsest lattice spacing is not included. The original and thermal-improved estimator are shown in the right-hand panel of Figure 6, which likewise strongly suggests the suppression of lattice artifacts in the difference. Once again, the difficulty of performing a controlled continuum limit purely based on the vacuum correlators at large virtualities is apparent even with fine lattice spacings available. For illustration, we quote an estimate for the improved observable where the central value and systematic error are determined in the same way as in the subsection before For the five-loop perturbative result, using and the renormalization scale µ = Q, we found D pert (Q 2 ) = 3.24(−3)(1). Clearly, the perturbative prediction is in good agreement with our non-perturbative estimate obtained with the help of thermal correlators computed at very fine lattice spacings.

V. SUMMARY AND OUTLOOK
We have shown that a better control over the continuum limit for short-distance dominated integrals over the vector correlator can be achieved by using thermal screening correlators at particularly fine lattice spacings, and then correcting for the difference. For the Adler function at a momentum Q = 2.36 GeV, we estimate that we achieved a reduction of the systematic error due to the continuum limit by a factor of four, relative to the continuum limit based on the available vacuum correlators only. This is significant, since taking the continuum limit is responsible for one of the leading systematic errors on this quantity. The cost of generating the thermal ensembles at small lattice spacings is modest in comparison to the cost of the vacuum ensembles. The method can also be applied to the charm contribution, which is even more short-range and therefore more susceptible to large cutoff effects.
For the anomalous magnetic moment of the muon, we find that, in the time-momentum representation, performing a 'naive' linear extrapolation of the t ≤ 0.2 fm contribution using lattice spacings down to 0.049 fm leads to an overestimate of about three percent as compared to our best estimate based on lattice spacings down to 0.033 fm. Since the muon (g − 2) is an observable which is much more infrared-weighted, the slightly inaccurate continuum extrapolation of the short-distance contribution in the (u, d, s) quark sector only amounts to a difference of about 0.4 × 10 −10 on this quantity, which is an order of magnitude smaller than the precision of the current most precise estimates [4,24]. Nevertheless, we have seen clear evidence that performing additively tree-level improvement reduces the lattice artifacts in this short-distance regime and we thus recommend its use in future calculations of the leading hadronic contribution to (g − 2) µ .
One may wonder, is it possible to calculate the hadronic contribution to the running of the QED coupling up to the Z-boson mass in lattice QCD with controlled errors? The methods and tests presented in this paper strongly suggest that it is indeed possible, and we now sketch a promising strategy. Let ∆ 2 Π(Q) ≡ Π(Q 2 ) − Π(Q 2 /4) be the difference of vacuum polarisations corresponding to the running of α over the momentum interval Q/2 to Q. Such an observable is very similar to the Adler function D(Q 2 ), in that it is dominated by Euclidean distances of order 1/Q. What we have seen above leads us to conclude that ∆ 2 Π(Q 2 ) for a large Q 1 GeV can be computed at the one-percent level at a temperature T ≈ Q/(8π) in terms of the screening vector correlator; we may write this quantity ∆ th 2 Π(Q; T ). The difference ∆ th 2 Π(Q; T /2) − ∆ th 2 Π(Q; T ) can be evaluated by performing simulations at temperatures differing by a factor of two with common lattice spacing. See Appendix E for an encouraging study at leading order in perturbation theory. Parametrically, the difference between ∆ 2 Π(Q) and ∆ th 2 Π(Q; T ) is of order (πT /Q) 4 , and thus ∆ th 2 Π(Q; T /2) − ∆ th 2 Π(Q; T ) already provides a sufficiently good estimate for that difference. The latter can also be estimated using perturbation theory, since a high value of Q is tied to a high value of T . Thus by varying the temperature by factors of two, we can map out the running of α by factors of two in Q up to the Z-boson mass. Such a program can be carried out on existing computing platforms, albeit at a significant investment, since realizing the double hierarchy πT Q π/a typically requires the use of lattices of size 48 × 192 3 . If one resorts to such fine lattices, it is probably mandatory to address the freezing of the topological charge, for instance by using open boundary conditions in the x 3 -direction [25,26].
better agreement between the vacuum and thermal correlation functions at short distance, due to the presence of constant terms which are not captured in the OPE picture. This point will be illustrated in the following.
As we discussed in some detail the OPE-based expansion of the electromagnetic-current correlator in a recent publication ( [27]), we largely refer to that for setting the notation and for the introduction of the main observables. More specifically, we make reference to the first few paragraphs of Section 3 for the definition of the Wilson coefficients and of the tower of local twist-two operators, and to the beginning of Section 3.1 for the operator mixing in the dimension-four sector. Moreover, we refer to Appendix B for the operator expectation values in the theory of free quarks. medium is at rest by setting u = (1, 0). Concentrating on the choices of indices µ = ν = 0, 1, we obtain Fourier-transforming with respect to the variable q 3 by using Eq. (B1) in Ref. [28], whose infrared regularization results in the q 3 -independent term being zero, we find for the coordinate-space version of Eq. (A6) Finally, to get to the isovector vector-current correlators in Euclidean space-time we drop the factor f Q 2 f and multiply by −1 the contribution to the correlator with two spatial indices. We obtain the OPE prediction to leading order in the gauge coupling In Appendix B we verify the correctness of the O(|x 3 |) term in the theory of free quarks. We observe that the term linear in |x 3 | cancels out in the sum G th 00 (x 3 ) + G th 11 (x 3 ). This fact is not sufficient to conclude that taking this combination improves the short-distance agreement with the vacuum correlation function G 00 (x 3 ) = G 11 (x 3 ), because the OPE is not able to capture constant terms in the small-x 3 expansion, as those are not of short-distance origin. The free-theory correlators provide a concrete example of this situation, as shown in Figure 7. The figure shows the difference between the thermal correlation functions G th 00 (x 3 ), G th 11 (x 3 ) and their vacuum counterparts as a function of x 3 , between x 3 = 0 and x 3 0.4/T . The leading OPE term, as explicitly computed in Appendix B, is also reported. Contrary to the OPE prediction, the difference between the thermal and vacuum correlators starts from a nonzero value in x 3 = 0. However, once we subtract ∆G 0 ≡ [G th − G] x 3 =0 , we find that the truncated OPE provides indeed a good description of (G th − G − ∆G 0 ) at small x 3 .

Next-to-leading-order Wilson coefficients
To higher order in perturbation theory, the mixing under renormalization with the gluonic twist-two operators O µ 1 ...µn ng (Eq. (3.3) in Ref. [27]) must be taken into account. To write explicitly the contribution from the dimension-four operators to NLO precision, we follow closely the sections 3 and 3.1 of Ref. [27] and refer to those for any unexplained notation. Going back to Eq. (A6), its NLO equivalent is obtained by making the substitution where e and p are the energy density and the pressure of the thermal medium andμ and Λ represent two energy scales, respectively the one at which the local operators are renormalized and the one at which the one-loop renormalized coupling diverges. The power γ is one of the eigenvalues of the anomalous-dimension matrix, where b 0 = 11 − 2 3 N f is the coefficient of the one-loop contribution to the beta function. To get the final expression in position space, we are faced with the problem of computing the Fourier transform 1 2π As in the case of Eq. (A7), the Fourier integral is divergent and a regularization procedure is needed to correctly extract the asymptotic dependence on x 3 . We discuss this point in detail in Appendix C and report here the final result The OPE prediction to NLO in the gauge coupling reads where µ = 3 and s µ is a sign factor evaluating to 1 for µ = 1, 2 and to −1 for µ = 0. While we did not discuss explicitly the case µ = 2 so far, we point out that the symmetries of the screening correlator constrain it to be equal to the case µ = 1. In the extreme x 3 → 0 limit the logarithmic contribution is subleading, and we find As already pointed out in Ref. [27] relative to the structure functions of the quark-gluon plasma, the free-theory prediction of Eqs. (A9), (A10) is not recovered in the short-distance limit x 3 → 0, contrary to the intuition coming from the property of asymptotic freedom. As a consequence of the mixing between operators under renormalization, the free theory represents here an extreme case which is not connected to the real interacting theory by an expansion in powers of the QCD coupling.
which is 7πN c T 4 90 1 r . (B6) To project to zero momentum in the directions x 1 and x 2 , we make use of the following equation To conclude, the linear term in the expansion of the correlator G th 00 (x 3 ) in powers of |x 3 | is given by and it is in agreement with the OPE prediction −2|x 3 |T 2 O 2f We project to zero momentum in the directions x 1 and x 2 by applying the following procedure To conclude, the linear term in the expansion of G th 11 (x 3 ) in powers of |x 3 | is 7π 2 N c T 4 45 and it agrees with the OPE prediction 2|x 3 |T 2 O 2f free T . Table VI: Accuracy of the continuum extrapolation of the thermal quantity ∆ th 2 Π(Q 2 ; T ), measured against the correct continuum value. The label "plain" refers to fitting the plain lattice observable, while for the case "subtr." the logarithmic lattice artifact is subtracted as ∆ th 2 Π(Q 2 ; T ) −c ∆ 2 Π a 2 log(1/(T a)). For the choice of scales made here T = 246.25 MeV, Q = 8πT 6.2 GeV, we havec ∆ 2 Π 2.2 fm −2 (see Eq. (E3)). The lattice data and the fit curves are shown in Figure 8, left panel.  6.2 GeV. The O(a 2 log(1/a)) lattice artifact is either included in the fit ansatz or explicitly subtracted by using the known prefactorc ∆ 2 Π (E3). The accuracy of the resulting continuum estimates is reported in Table VI. Right: continuum extrapolation of the bias correction [∆ th 2 Π(Q 2 ; T /2) − ∆ th 2 Π(Q 2 ; T )], together with its known continuum value. The fit ansatz is linear in a 2 .
To conclude, the free-theory analysis indicates that the strategy outlined in Section V to compute the hadronic contribution to the running of α em up to large energy scales (compared to the hadronic scales easily accessible on the lattice) is feasible with the setup and lattice sizes suggested there. The idea is to study the difference of hadronic vacuum polarizations ∆ 2 Π(Q 2 ) between the scales Q and Q/2 making use of thermal lattices with temperatures T ≈ Q/(8π) and T /2. Having fixed T = 246.25 MeV and Q = 8πT 6.2 GeV, and having considered lattices with at the most 48 points in the Euclidean-time direction, we