The hadronic vacuum polarization contribution to the muon $g-2$ from lattice QCD

We present a calculation of the hadronic vacuum polarization contribution to the muon anomalous magnetic moment, $a_\mu^{\mathrm hvp}$, in lattice QCD employing dynamical up and down quarks. We focus on controlling the infrared regime of the vacuum polarization function. To this end we employ several complementary approaches, including Pad\'e fits, time moments and the time-momentum representation. We correct our results for finite-volume effects by combining the Gounaris-Sakurai parameterization of the timelike pion form factor with the L\"uscher formalism. On a subset of our ensembles we have derived an upper bound on the magnitude of quark-disconnected diagrams and found that they decrease the estimate for $a_\mu^{\mathrm hvp}$ by at most 2%. Our final result is $a_\mu^{\mathrm hvp}=(654\pm32\,{}^{+21}_{-23})\cdot 10^{-10}$, where the first error is statistical, and the second denotes the combined systematic uncertainty. Based on our findings we discuss the prospects for determining $a_\mu^{\mathrm hvp}$ with sub-percent precision.


Introduction
After the discovery of the Higgs boson the search for physics beyond the Standard Model has further intensified. The three principal strategies include the observation of new particles, the detection of enhanced signals in rare decay processes and deviations between experimental determinations of precision observables and theoretical predictions based on the Standard Model. One of the most prominent examples for the latter is the value of the anomalous magnetic moment of the muon, a µ = 1 2 (g − 2) µ , which exhibits a persistent deviation of 3.6σ at the current precision of 0.5 ppm [1]. It is well known that the theoretical uncertainty is dominated by hadronic contributions, more specifically the hadronic vacuum polarization and hadronic light-by-light scattering contributions, a hvp µ and a hlbl µ , respectively. The estimate for a hvp µ which enters the Standard Model prediction is typically obtained from dispersion theory using the experimentally determined cross section e + e − → hadrons as input [2][3][4][5][6][7].
In this paper we present results for a hvp µ in lattice QCD, using two complementary approaches: The first is based on the standard determination of the vacuum polarization function Π(Q 2 ) via a four-dimensional Fourier transform of the vector correlator. The second method uses the so-called "time-momentum representation" (TMR) discussed in [51,54,57]. As another variant we consider time moments of the vector correlator [34] to describe the low-momentum region of Π(Q 2 ). We focus primarily on controlling the various sources of systematic uncertainties associated with the lattice approach to a hvp µ , and in particular the problem of constraining the deep infrared region.
Our work is based on QCD with two light degenerate dynamical quarks. The inclusion of the effects from isospin breaking and from dynamical s, c and b quarks is left for future work. Clearly, for a precision determination of a hvp µ in lattice QCD it is necessary to include dynamical strange and charm quarks. However, the collection of results for a wide range of quantities in [58] suggests that the effects from the strange and charm quarks in the sea can be expected to be subleading at our level of precision. While the calculation of quarkdisconnected diagrams has only been performed on a subset of our ensembles, this has still allowed us to derive an upper bound on their overall influence which is included in the final error estimate. Our main result, stated in eq. (37), is the determination of a hvp µ with an overall precision of 6%. While this is still significantly larger than the quoted uncertainty of the dispersive approach, our study provides valuable insights for future lattice calculations of this important quantity.
This paper is organized as follows: In section 2 we discuss different approaches for computing the hadronic vacuum polarization contribution to (g − 2) µ . Simulation details are described in section 3, and in section 4 we present a detailed discussion and comparison of our results obtained on individual ensembles. The extrapolation of our results to the physical point is described in section 5, including a detailed discussion of systematic errors. We state our conclusions in section 6. In a series of appendices we present further details on the current renormalization, the efficient evaluation of the QED kernel in the TMR, the estimation of finite-volume effects and the calculation of quark-disconnected diagrams, respectively.
2 Lattice approaches to a hvp µ The hadronic vacuum polarization contribution, a hvp µ , to the muon anomalous magnetic moment can be obtained from the vacuum polarization function Π(Q 2 ) convoluted with a known kernel function K(Q 2 ; m 2 µ ) (defined in appendix B) and integrated over Euclidean momenta Q 2 [28,50,59], as where α and m µ are the electromagnetic coupling and muon mass, respectively. The vacuum polarization function Π(Q 2 ) is obtained from the vacuum polarization tensor Π µν (Q), which is given in terms of the correlator of the electromagnetic current J µ (x) as where Q denotes the Euclidean momentum. Euclidean O(4) invariance and current conserva- The subtracted vacuum polarizationΠ(Q 2 ), defined bŷ which appears in the integrand, is free of UV divergences. Using the explicit expression for the kernel function [28,60] one infers that the integrand in eq. (1) is peaked near Q 2 ≈ m 2 µ ≈ 0.01 GeV 2 . To access such small momenta on a finite lattice directly would require volumes corresponding to a linear extent of O(10 fm) or more, which is difficult to achieve with currently available resources. Therefore, the exact shape of Π(Q 2 ) in the small-momentum region, as well as the value of Π(0) are difficult to determine with sufficient accuracy.
Several methods for accurately constraining the small-momentum regime have been proposed and studied. This includes the use of twisted boundary conditions [61][62][63] that are designed to penetrate more deeply into the region near Q 2 = 0 [32,64,65], and the direct determination of the additive renormalization Π(0), either via operator insertions [53] or by computing time moments of the vector correlator [34]. In order to avoid introducing any model dependence it has been proposed to represent Π(Q 2 ) by either Padé approximants or conformal polynomials in a sub-interval 0 ≤ Q 2 ≤ Q 2 cut and to evaluate the convolution integral for momenta Q 2 > Q 2 cut using the trapezoidal rule [56]. Such a "hybrid strategy" requires accurate data for sufficiently small values of Q 2 cut . In the so-called "time-momentum representation" (TMR) discussed in [51,54,57] the subtracted vacuum polarization functionΠ(Q 2 ) is directly obtained from the spatially summed two-point correlator G(x 0 ) of the electromagnetic current, i.e.
When inserted into eq. (1), the hadronic vacuum polarization a hvp µ is given by where the x 0 -dependent kernel function K(x 0 ; m µ ) is obtained by performing the integral and K(Q 2 ; m 2 µ ) is the same kernel function as in eq. (1). A representation of K(x 0 ; m µ ) suitable for a numerical evaluation is given in appendix B. The main technical difficulty in this approach arises from the fact that the vector correlator G(x 0 ) is integrated to infinite Euclidean time. Therefore, the large-x 0 behaviour of G(x 0 ) must be accurately constrained. For light enough pion masses the vector correlator is dominated by the two-pion state as x 0 → ∞, and thus one has to resort to elaborate calculations of G(x 0 ) including multi-particle states [51].
A closely related method for determining the subtracted vacuum polarization function Π(Q 2 ) is based on the calculation of the time moments of the vector correlator [34]. The starting point is the expansion of Π(Q 2 ) at low Q 2 , i.e.
When Q is chosen as Q = (ω, 0) one obtains the vacuum polarization function (VPF) from the spatially summed vector correlator G(x 0 ) according to The expansion coefficients Π 0 , Π 1 , Π 2 , . . . in eq. (9) can be determined from the derivatives with respect to ω 2 which are, in turn, related to the time moments G 2j of the vector correlator via In this way one obtains The time moments can be used to construct the Padé representation of the subtracted VPF Π(Q 2 ) ≡ 4π 2 (Π(Q 2 ) − Π(0)) in the low-momentum regime. There is also a close relation between time moments and the TMR: by expanding the sine function in eq. (5) as a power series in Q 2 one recovers the time moments as expansion coefficients in accordance with eq. (9). For later use it is also convenient to consider the decomposition of the electromagnetic current into an iso-vector (I = 1) and an iso-scalar (I = 0) part, according to where we use the superscript ρ to denote the iso-vector (I = 1) contribution. The corresponding correlator is defined by and the iso-spin decomposition of the vector correlator reads Note that only quark-connected diagrams contribute to the iso-vector correlator G ρρ (x 0 ).

Simulation details
Our calculations have been performed on a set of ensembles with N f = 2 flavours of dynamical, mass-degenerate, O(a)-improved Wilson quarks and the Wilson plaquette action. The improvement coefficient c sw was tuned according to the non-perturbative determination of ref. [66]. The gauge configurations have been generated as part of the CLS (Coordinated Lattice Simulations) initiative, using the deflation-accelerated DD-HMC [67,68] and MP-HMC [69] algorithms.
In Table 1 we have compiled the parameter values, system sizes and overall statistics used in our determination of the hadronic vacuum polarization contribution. The values for the lattice scale reported in the table have been determined using the kaon decay constant [70,71]. In order to enhance statistics we have used four source positions per configuration, except for the most chiral ensembles G8 and O7 for which up to 16   resulting number of measurements for each ensemble is shown in the right-most column of Table 1.
The bare values of the strange quark mass used in this work are based on an update of the analysis of ref. [70] where the physical values of the kaon mass and decay constant were used to set κ s . The updated analysis [72] includes improved determinations of the renormalization factors Z A of the axial current, increased statistics, as well as a new measurement of κ s for the ensembles B6 and G8. In the charm sector, we used the bare quark masses determined from the experimental value of the D s -meson mass in ref. [73] for the two finest values of the lattice spacing. Based on these results, at β = 5.2 we estimated the hopping parameter κ c of the charm quark from the a 2 dependence of the ratio, m c /m s . Values for κ s and κ c are listed in Table 2.
In our calculation we have considered a mixed vector correlator including the conserved point-split vector current and the local vector current where f denotes one of the quark flavours u, d, s and c. The local current is neither conserved nor improved, yet it can be renormalized in a fashion that is consistent with O(a) improvement Here m f denotes the bare subtracted quark mass of quark flavour f , b V and c V are improvement coefficients, and T µν, is the tensor current. The conserved vector current, while not subject to renormalization, requires O(a) improvement even at tree level, which was not considered in this work. Since we did not determine the matrix elements containing the derivative of the tensor current, our results for a hvp µ are not fully O(a) improved.
In the light quark sector the mass-dependent factor in eq. (18) is usually a small correction. However, since we compute the contribution from the charm quark to a hvp µ , the corresponding mass dependence is sizeable and must be included for a reliable extrapolation to the continuum limit. We have considered two different procedures for the determination of the renormalization factor of the local vector current, including the mass dependence: 1. Determine Z V using the interpolating formula in ref. [75] and evaluate the one-loop expression for the improvement coefficient b V from [76] using the boosted coupling g 2 ≡ g 2 0 / 1 3 Tr U P . as determined via the second procedure. As will be discussed in detail in section 5, we observe large lattice artefacts in the case of the charm quark contribution to a hvp µ . In order to check for the stability of the continuum extrapolation we have compared the results obtained using both procedures to determine the current normalization and found very good agreement.
With the above definitions of the currents, the vacuum polarization tensor can be expressed in terms of the mixed vector correlator as where q f , q f denote the electric charges of quark flavours f and f , andQ µ = 2 a sin aQµ 2 is the lattice momentum. Like in our previous publication [32] we have used twisted boundary conditions [61][62][63] in order to apply an additive shift to the momentum of the quark propagator. In this work we used a single value of the twist angle, chosen such as to provide three equidistant values of Q 2 between the lowest two Fourier momenta, as well as one additional data point below (2π/L) 2 . The imposition of twisted boundary conditions induces the breaking of isospin symmetry and modifies the Ward identity of the vacuum polarization tensor that guarantees its transversality [64]. We have checked explicitly [77] that the violation of the Ward identity has a negligible effect on the determination of Π(Q 2 ).
It has been noted in [51,78] (see also [44,79]) that the vacuum polarization tensor does not vanish at Q = 0 in finite volume, Π µν (0) = 0. In order to reduce finite-volume effects it is then advantageous to subtract the contribution Π µν (0), which is easily effected via a simple modification of the phase factor in eq. (19), i.e.
Run am π am ρ m π /m ρ κ s am V (ss) κ c am V (cc) A3 0.1893(6) 0.3937 (29)   In addition to computing Π µν (Q) we have also considered the spatially summed vector correlator, given by The sum f,f . . . in equations (19) and (21) runs over all quark flavours included in the electromagnetic currents. Here we focus on the quark-connected contributions to the vector correlator. In order to quantify individual flavour contributions to a hvp µ it is useful to define where q 2 ud = 5/9, and it is understood that the expectation value is restricted to quarkconnected diagrams. The vector correlator in the long-distance regime is constrained by the mass spectrum of the theory. Depending on the value of the light quark mass on a given ensemble, the lowest-lying state corresponds either to the vector meson or to a twopion state. For a reliable determination of the energy levels in the vector channel, we have computed additional correlators using standard Gaussian smearing [80] in the calculation of quark propagators, with APE-smeared link variables [81] in the spatial directions. The mass in the vector channel and also the pion mass used in this study were determined from the appropriate correlation functions with smearing applied both at the source and sink. The corresponding mass estimates are listed in Table 2.
All statistical errors were estimated using a bootstrap procedure with 10,000 samples. For the estimation of systematic errors we employed the so-called "extended frequentist method" [82,83] and determined the distributions of results obtained from a set of variations of our analysis procedure. Details are provided in the sections describing our results.

Calculation of a hvp µ
In this section we report on the determination of a hvp µ for all our ensembles, employing different methods, in order to check for systematic effects.

a hvp
µ from the hybrid method Our calculation of a hvp µ from the vacuum polarization tensor proceeds by evaluating the vacuum-subtracted tensor defined in eq. (20) and factoring out the tensor structure according to eq. (3). In order to determine the additive renormalization Π(0) and describe the data in the small momentum regime, we have employed the ansatz where P [n,m] denotes the Padé approximant of order [n, m]. Following ref. [52] we consider n = m or n = m + 1 and write P [n,m] as In accordance with the discussion of the "hybrid strategy" in [56] the main task is to determine the Padé representation in an interval 0 < Q 2 < ∼ Q 2 cut . Here we have adopted two procedures: the first proceeds by determining the coefficients A k and B k from fits to the VPF, the second uses time moments to construct the Padé approximation for 0 < Q 2 < ∼ Q 2 cut . Ideally, the Padé representation of Π(Q 2 ) should be constructed by considering a sequence of approximants of increasing order [52]. However, when confronted with actual simulation data one often finds that the latter are not constraining enough to allow for a systematic investigation whether successive Padés converge towards the actual VPF. One therefore resorts to constructing low-order Padé approximations, i.e. one-pole ansätze that are not much different from a vector meson dominance description. To minimize the bias incurred from using a particular Padé approximant, the value of Q 2 cut should be chosen much smaller than m 2 ρ . However, one has to balance this requirement against fit stability and statistical accuracy. In order to have sufficiently many data points available so that stable correlated fits with acceptable χ 2 /dof can be performed, we have chosen Q 2 cut ≈ 0.5 GeV 2 . At our level of statistical precision we find that the data are well described by a Padé [1,1] ansatz and exhibit values of the correlated χ 2 /dof of order unity, except for ensembles E5 and N6 for which χ 2 /dof > 6. Using a Padé [2,1] ansatz gave consistent results but larger statistical errors.
In order to calculate the light quark contribution to the anomalous magnetic moment, (a hvp µ ) ud , we have evaluated the convolution integral of eq. (1) in the interval 0 ≤ Q 2 ≤ Q 2 cut by inserting Π(Q 2 ) ud − Π(0) ud as determined by the Padé [1,1] fit. The contribution from the region Q 2 > Q 2 cut was computed using trapezoidal integration, and the resulting values of (a hvp µ ) ud are shown in the third column of Table 3. To check for stability against variation of the scale Q cut we have computed (a hvp µ ) ud for Q 2 cut ≈ 0.3 − 0.35 GeV 2 . We find agreement within slightly larger errors with the numbers reported in Table 3.
For the determination of the strange quark contribution to the vacuum polarization, Π(Q 2 ) s − Π(0) s , and the anomalous magnetic moment, (a hvp µ ) s , we have followed the same  procedures as for (a hvp µ ) ud . Concerning the influence of variations in the value of Q 2 cut and the order of the Padé ansatz we came to the same conclusions. The results for (a hvp µ ) s determined for Q 2 cut ≈ 0.5 GeV 2 are listed in the fourth column of Table 3. For ensemble E5 we again found χ 2 /dof ≈ 7, both for the Padé [1,1] and [2,1] fits. The corresponding entry is marked by an asterisk in Table 3 and is excluded from the subsequent analysis.
The Q 2 -dependence of the charm quark contribution to Π(Q 2 ) shows a lot less curvature compared to the lighter flavours. We have therefore applied a slightly different procedure, by fitting Π(Q 2 ) not only to a Padé [1,1] ansatz but also to a linear function in Q 2 . Starting from Q 2 cut ≈ 0.5 GeV 2 we have gradually lowered Q 2 cut until the two different ansätze gave consistent results. The corresponding estimates of (a hvp µ ) c are listed alongside with the respective values of Q 2 cut in Table 3. A striking but not unexpected feature of (a hvp µ ) c is the strong dependence on the lattice spacing. This is seen easily by comparing the estimates for (a hvp µ ) c for ensembles B6, F7 and O7: at approximately constant pion mass in physical units the results for (a hvp µ ) c vary by 30-40% within the range of lattice spacings considered in this work.
An alternative determination of the low-momentum representation ofΠ(Q 2 ) is achieved by computing time moments of the vector correlator. These are linked to the coefficients Π j in the Taylor-series expansion of the vacuum polarization function and also to the additive renormalization Π(0) (see eq. (12)). The Π j 's can then be used to construct the coefficients A k , B k in the Padé representation of eq. (24). For instance, the Padé [1,1] approximant written in terms of the expansion coefficients reads and expressions for higher-order Padés can easily be worked out. The determination of the time moments proceeds by summing the vector correlator over all Euclidean times. As in the case of the TMR, which is discussed in detail in the next subsection, this requires some sort of modelling of the long-distance regime of G f (x 0 ). To this end we have assumed that G f (x 0 ) is described by a single exponential for x 0 > x cut 0 (see eq. (28) below). A more detailed discussion is presented in section 4.2.
It is instructive to compare the Padé representation of Π(Q 2 ) as determined from time moments to that obtained from fits to Π(Q 2 ) below Q 2 cut discussed earlier. Such a comparison is shown in Fig. 1 for the ensembles G8 and O7. In particular, we compare the intercept Π(0) ud as obtained from a Padé [1,1] fit for 0 < Q 2 ≤ 0.5 GeV 2 to its determination from the second time moment. As is apparent from the figure the two procedures agree very well, which is an important cross check. Typically, the estimate of Π(0) ud from the fit has a smaller error. Having computed the coefficients Π 0 , Π 1 , . . . , Π 4 from time moments we constructed the Padé [1,1] and [2,1] representations ofΠ(Q 2 ) in the interval 0 ≤ Q 2 ≤ Q 2 cut . As before we determined a hvp µ by performing the convolution integral overΠ(Q 2 ) for Q 2 > Q 2 cut using trapezoidal integration. Thus, our way of employing time moments differs from the procedures applied in refs. [34,39], where the subtracted vacuum polarization functionΠ(Q 2 ) is constructed from time moments within the entire momentum interval.
In order to guarantee a smooth transition between the low-momentum representation and the actual data forΠ(Q 2 ) = 4π 2 (Π(Q 2 ) − Π 0 ) we have chosen Q 2 cut so as to minimize  the difference between the Padé approximation ofΠ(Q 2 ) and the data within the interval Q 2 = 0.1 − 0.5 GeV 2 . Results for a hvp µ obtained via this procedure are listed in Table 4. We found the differences between the Padé [1,1] and [2,1] descriptions of the low-Q 2 regime to be negligible.

The TMR method for a hvp µ
The integral representation of the subtracted vacuum polarization function,Π(Q 2 ), is shown in eq. (5), and the hadronic vacuum polarization contribution of quark flavour f = (ud), s, c to a µ is then obtained as [51], In appendix B we derive an explicit expression which describes K(x 0 , m µ ) with an accuracy of O(10 −6 ). The kernel is proportional to x 4 0 at small x 0 , and to x 2 0 at large x 0 . The integration must be performed over all Euclidean times x 0 , and thus the challenge in this method is to control the long-distance behaviour of the spatially summed vector correlator G f (x 0 ) defined in eq. (22). The main issues are that In order to handle the large-x 0 part separately, we define our estimator The subscript "inter" denotes that the vector correlator has been obtained from a local cubic spline interpolation of the numerical data. The long-distance part G f (x 0 ) ext is obtained by extending the correlator by one of the methods specified below. Items (a) and (b) can be dealt with by extrapolating the correlator using a sum of exponentials. Indeed, in a finite volume, the spectral representation implies that the correlator is exactly given by an infinite sum of exponentials exp(−E n x 0 ). The lowest few energyeigenstates 1 dominate at large x 0 . Therefore the simplest incarnation of this method is to use a single-exponential extension of the correlator, The parameters (A, m V ) depend on the flavour composition f = (ud), s, c of the vector current. Clearly, the systematic error incurred by using a single exponential must be investigated. Since the energy levels only depend on the quantum numbers of the interpolating operator, they can also be determined from auxiliary correlation functions. In our benchmark analysis, whose preliminary results have been presented in [84], we extract m V from the two-point function of a smeared vector operator, obtaining the masses reported in Table 2. The amplitudes A are then determined from a one-parameter fit to eq. (28) using these masses as input. A compilation of results for a hvp µ extracted via the TMR is shown in Table 5 along with the respective values of x cut 0 . As an illustration of the method, we plot the integrand of eq. (26) for the light-quark connected contribution on the two ensembles with the lightest pion masses, G8 and O7, in Fig. 2. The extension method just described is labelled as '1-exp'. Various coloured bands represent other methods (discussed below) to constrain the long-distance behaviour of the vector correlator.
The choice of x cut 0 affects the accuracy of a hvp µ since larger values of x cut 0 increase the statistical error because of the quickly rising noise-to-signal ratio in the correlator data. By contrast, a smaller cutoff implies that estimates of a hvp µ will be more strongly affected by systematic effects arising from assumptions regarding the asymptotic behaviour of the correlator. We have chosen x cut 0 as the value beyond which the statistical signal deteriorates to such an extent that the original data do not accurately constrain the correlator anymore. In terms of statistical accuracy this represents the most conservative choice, since the interpolation of G f (x 0 ) is used within the maximum Euclidean time range where the signal is not lost. We have checked explicitly that our estimates are not affected by the particular choice of x cut 0 . Moreover, in the case of the strange and charm quark contributions we have found that the correlators fall off so rapidly that the effect of truncating the integral in eq. (26) at x 0 = x cut 0 on the estimates of (a hvp µ ) s and (a hvp µ ) c is insignificant. We conclude that in this case the systematic error arising from the modelling of the long-distance contribution is negligible for x cut 0 > ∼ 1.2 fm. In the future, variance-reduction strategies, such as those described in [85,86] may be used to suppress the strong growth of the noise-to-signal ratio of G f (x 0 ), thereby reducing the need for modelling the large-x 0 behaviour.
We now return to the issue of the extension of the correlator G ud (x 0 ). On all our ensembles except for G8, a single exponential already provides a remarkably good description of the correlator for x 0 ≥ x cut 0 . The reason is that the lightest energy-eigenstate in the box has a  that marks the switch from a cubic spline interpolation of the correlator to its long-distance representation. The label "1-exp" refers to the single exponential of eq. (28), while "GS" and "GS, inf" refer to the Gounaris-Sakurai-based extensions in finite and infinite volume, respectively. For the latter a slightly smaller value of x cut 0 was used on ensemble G8 to stabilize the fit. At heavy pion mass only the one-exponential extension was considered.
large amplitude relative to the other states. This fact is well understood: The finite-volume energies and amplitudes are directly related to the timelike pion form factor [87,88]. The latter peaks at the ρ-resonance, E = m ρ , and one state in the finite box almost always lies nearby in energy. It happens to be the lightest state on all but one ensemble. Thus the reason that the light-quark correlator G ud (x 0 ) is dominated by a single exponential is closely related to the ideas underlying the vector-meson dominance model (VMD) used in hadron phenomenology.
Obviously, the one-exponential extension has its limitations. This becomes most evident on ensemble G8, where one expects to find, below the energy level E 2 associated with a large amplitude, an energy level E 1 < E 2 with a smaller amplitude. This conclusion is easily reached by initially neglecting the interactions between two pions in the T 1 representation, E ππ ≡ E 1 = 2 m 2 π + (2π/L) 2 (≈ 695 MeV on G8). The non-vanishing scattering phase leads to a modest shift of the energy level. Obviously the result for a hvp µ incurs a bias if one ignores this low-lying state, but it is difficult to determine its precise energy and amplitude from G ud (x 0 ), because the amplitude is small. These observations also show that the finitevolume correlator behaves drastically differently at large x 0 than in infinite volume: in the latter case, G ud (x 0 ) is dominated by a two-pion continuum starting at E = 2m π (≈ 370 MeV on G8) rather than by discrete energy levels. Thus the issue of extending the correlator G ud (x 0 ) to long distances is intimately related to the question of the finite-size effects on lattice determinations of a hvp µ (see item (c) above). To prepare for a more sophisticated treatment of the long-distance behaviour of the vector correlator, it is useful to recall the isospin decomposition of eq. (15), i.e. G(x 0 ) = G ρρ (x 0 ) + G I=0 (x 0 ). The iso-vector part G ρρ is directly proportional to the quark-connected light-quark contribution G ud , i.e.
The ω-resonance is the lowest-lying state in the iso-scalar channel, which has a much smaller width compared to the ρ. In particular, the decay of the ω into three pions is strongly suppressed, and thus the single exponential is a good approximation for evaluating the iso-scalar contribution to the convolution integral in eq. (7). By exploiting the fact that the ρ − ω splitting is small, we arrive at our final ansatz for the long-distance contribution to the quark-connected light quark vector correlator, i.e.
In other words, we replace the light iso-scalar correlator by a single exponential with m V = m ρ in the long-distance regime. 2 In the following subsection we describe how G ρρ (x 0 ) ext can be constrained via the Gounaris-Sakurai model.

Gounaris-Sakurai based extension of the vector correlator
As already advocated in [51], the calculation of the vector correlator for a hvp µ should ideally be accompanied by a dedicated study of the timelike pion form factor F π (ω). This has been the subject of a few recent lattice calculations [89][90][91]. With the pion form factor at hand, the long-distance part of the iso-vector correlator G ρρ (x 0 ) ext can be obtained straightforwardly. Moreover, one can compute the infinite-volume iso-vector correlator via thus correcting model-independently for the dominant finite-size effects in a hvp µ . Eq. (32) assumes that the 2π channel saturates the iso-vector correlator, which is a good approximation if x cut 0 is sufficiently large. However, lacking a full-scale calculation of the timelike pion form factor, we apply a simplified version of this strategy (at the cost of a certain modeldependence). Based on the success of the Gounaris-Sakurai (GS) model [92] in describing experimental data for e + e − → π + π − data, we assume that the timelike pion form factor is well approximated by this model at the pion masses used in our ensembles. Since the GS model only contains two parameters (the ρ-mass and its width Γ ρ ), the same number as the one-exponential ansatz eq. (28), this simplified approach allows us to go beyond the one-exponential extension whilst remaining numerically viable given the available lattice data. The procedure can be summarized as follows: 1. Fix the GS parameter m ρ by identifying it with one of the energy levels determined from the smeared-smeared correlator.
2. Determine the GS parameter Γ ρ from the iso-vector correlator G ρρ (x 0 ), using m ρ as input.
3. Determine the low-lying energy levels and their amplitudes using the GS model and the Lüscher formalism. The finite-volume correlator can then be computed beyond x cut 0 as the sum of the corresponding exponentials, and from there (a hvp µ ) ud is obtained.
4. In addition, the correlator G ρρ (x 0 ) can be calculated in infinite volume beyond x cut 0 via eq. (32), and from there a hvp µ is obtained. This estimator corrects for the dominant finite-size effects.
A discussion of the systematic error associated with the procedure is presented in appendix C. In steps 3 and 4, the lattice data G f (x 0 ) inter is used directly up to x cut 0 . We remark that the parameters describing the pion form factor must, in general, be determined simultaneously from the spectrum and finite-volume matrix elements; however in the present case we exploit the fact that the lowest two energy levels are only weakly dependent on Γ ρ .
To determine the GS ρ-mass from the smeared-smeared correlator (step 1) we have proceeded in the following way. For ensembles O7, N6, F7, F6, B6 and A5, the ρ-mass parameter of the GS model was extracted from a single-exponential fit to the smeared-smeared correlator. We have checked in these cases that, if the form factor is described by the GS model, identifying the lowest-lying energy-level with the GS ρ-mass is an excellent approximation, almost irrespectively of the value of Γ ρ . On ensemble G8, we have applied a two-exponential fit where the first energy level is set to E 1 = 2 m 2 π + (2π/L) 2 by hand and the second exponential is fitted and its mass identified with m ρ . In addition, both amplitudes A 1 and A 2 are fitted. Even with this three-parameter fit, we encountered a few bootstrap samples where the fit was unstable. Therefore we stabilized the fit in the following way: based on the ensembles with m π < 400 MeV, we performed an extrapolation of the GS ρ-mass linearly in m 2 π to the pion mass of the G8 ensemble, resulting in m xtrap ρ = (797 ± 15) MeV. We then used this information as a Bayesian prior, adding ∆χ 2 = (m ρ − m xtrap ρ ) 2 /σ 2 to the χ 2 , where σ was varied between 15 and 120 MeV. We found that the fit result was stable as long as σ ≤ 60 MeV. Fig. 2 shows the effect of describing the long-distance part of the correlator using the GS model as compared to using a single exponential for ensembles O7 and G8. Both the finitevolume and the infinite-volume versions are displayed. While the differences do not seem very dramatic, their impact on a hvp µ is significant, particularly because the effect of the two-pion continuum increases as the chiral limit is approached. By inserting the GS-based extensions of the iso-vector correlator G ρρ for x 0 > x cut 0 in finite and infinite volume into eq. (31) one can compute the corresponding estimates of (a hvp µ ) ud . The results are summarized in table 5.

Comparison of a hvp µ
We are now in a position to compare the estimates for a hvp µ obtained from different procedures described in the previous subsections. Obviously, this comparison refers only to the data without finite-volume corrections, since the latter have only been quantified for the TMR. The results listed in Tables 3, 4 and 5 show certain trends regarding their statistical errors. For instance, all three methods yield comparable statistical accuracy for the strange quark contribution (a hvp µ ) s . The light quark contribution (a hvp µ ) ud is equally precise when determined via the TMR or via Padé [1,1] fits below Q 2 cut . By contrast, constraining the low-Q 2 behaviour via time moments yields much smaller errors for (a hvp µ ) ud . Finally, the TMR is statistically by far the most precise method for determining the charm quark contribution (a hvp µ ) c . One might expect the results obtained using all three variants to agree for each individual ensemble. However, it is easy to see from Tables 3-5 that this is not always the case. The largest differences, which amount to about 10%, are observed for the charm quark contribution. By contrast, one mostly finds agreement among the estimates for (a hvp µ ) ud at the level of one or two standard deviations. Another interesting observation is the fact that the differences among estimates determined via the three methods decrease at smaller lattice spacing. Thus, the spread of results among individual ensembles can be attributed to a large part to the presence of lattice artefacts. This interpretation is further supported by the observationdiscussed in the next section -that the estimates for a hvp µ at the physical point agree within the quoted uncertainties.

Chiral and continuum extrapolations
We now describe our procedure for determining a hvp µ at the physical point, i.e. for vanishing lattice spacing and at the physical pion mass. We start by noting that there is no theoretically preferred ansatz which describes the chiral behaviour of a hvp µ in the range of pion masses which is usually considered in lattice simulations. We have therefore subjected the sets of results listed in Tables 3, 4 and 5 to simultaneous chiral and continuum extrapolations, using a variety of functional forms that parameterize the dependence on the pion mass and the lattice spacing, i.e. Fit A: Fit B: Fit C: Fit D: with fit parameters α 1 , α 2 , . . . , δ 2 . All four ansätze contain a term of order a, since the operators whose matrix elements determine the vacuum polarization are not fully O(a) improved.
The terms proportional to m 2 π ln m 2 π and m 4 π in fits A and B, respectively, account for the curvature in the chiral behaviour of the light-quark contribution (a hvp µ ) ud . By contrast, the pion mass dependence of (a hvp µ ) s and (a hvp µ ) c is mostly linear or even constant, which motivates the absence of such terms in fits C and D.
In order to estimate systematic errors associated with variations of our fitting and analysis procedures we have employed the so-called "extended frequentist's method" (EFM) [82,83]. When combined with the bootstrap method designed for the estimation of statistical errors one obtains the fit result from the median of the joint distribution, while statistical and systematic errors are represented by the lower and upper bounds of the central 68%. An overview of all fitting and analysis variants which enter the EFM are presented in Table 6. As regards variations of the ansatz for the chiral fit, we note that two additional functional forms were discussed in ref. [93], namely a fit including one inverse power of m 2 π , as well as a ChPT-inspired function containing a term proportional to ln m 2 π (i.e. without the factor of m 2 π multiplying the logarithm). We note that an ansatz containing ln m 2 π has a compelling justification only for m π < m µ [93] and does not apply to the situation realized in our simulations. While an inverse power of m 2 π does arise in the slope of Π(Q 2 ) at Q 2 = 0 via the numerically subdominant pion loop contribution [50], it may over-amplify the dependence of a hvp µ on m 2 π near the physical pion mass [93]. We have therefore excluded terms like 1/m 2 π and ln m 2 π from our EFM analysis. As a further check we have performed tentative fits based on a modified version of fit A, in which α 3 m 2 π ln m 2 π was replaced by α 3 ln m 2 π . The resulting estimates for a hvp µ at the physical point are well within the total error obtained by the EFM procedure. Thus, we conclude that the uncertainty associated with the chiral extrapolation has been quantified reliably.
The systematics of the chiral and continuum extrapolation can be investigated by varying the fit ansatz and by imposing different cuts in the maximum pion mass and the lattice spacing a. Another important systematic effect is associated with constraining the deep infrared regime of the vacuum polarization: In the case of the hybrid method we have used different TMR light strange charm Fit ansatz A, B A, B, C C, D Cuts in m π no cuts and a cut 1 cut 2 cut 2 cut 2 cuts 1 and 2 cuts 1 and 2 cuts 1 and 2 IR regime single exponential ‡ single exponential single exponential MeV † cut 2: a < 0.07 fm ‡ single exponential is not used as a variation with the GS model including the FV correction Table 6: Overview of variants of the fitting and analysis procedures which enter the estimation of systematic errors via the extended frequentist method. We focus on the hybrid method with the low-Q 2 behaviour determined by fits, as well as the TMR. The meaning of the various cuts is explained below the table. values of the momentum scale Q 2 cut below which the vacuum polarization function is described by a low-order Padé approximant.
For the TMR we have included two different variants for extending the vector correlator G ud (x 0 ) beyond x cut 0 , the first being the single-exponential ansatz, with the GS model (excluding the finite-volume correction) as an alternative. The GS-parameterization including the finite-volume shift was extrapolated separately. In this case we did not study effects of another ansatz for describing the infrared behaviour. For the strange and charm quark contributions we only used the single-exponential extension, since the estimates for (a hvp µ ) s and (a hvp µ ) c do not depend strongly on the details of the corresponding vector correlators for x 0 > ∼ 1.2 fm. The contribution from the charm quark to a hvp µ is particularly sensitive to the discretization and renormalization effects. This can be inferred already from the fact that the estimates for (a hvp µ ) c differ by 30-40% between our coarsest and finest lattice spacing (see Tables 3-5). Furthermore, combined chiral and continuum fits of the data including all three lattice spacings produce large values of χ 2 /dof, which is particularly pronounced for the data obtained  for the hybrid (above) and TMR (below) methods. Yellow bands correspond to the chiral behaviour in the continuum limit, while the dark red and blue curves represent the pion mass dependence at β = 5.5 and 5.3. The physical value of the pion mass is indicated by the vertical lines.
using the TMR. We have therefore consistently excluded the TMR-data for (a hvp µ ) c computed at the coarsest lattice spacing from the extrapolations to the physical point. Furthermore, in order to study whether the details of fixing the renormalization factor of the local vector current have a noticeable systematic effect on the extrapolation we have repeated the fits of (a hvp µ ) c using the factor . Another comment on the use of time moments to constrain the low-Q 2 dependence of Π(Q 2 ) is in order. We found that the combined fits to the results listed in Table 4 produced values of χ 2 /dof between 5 and 10 , regardless of the fit ansatz or of any other procedural variation. The most likely explanation is the smallness of the statistical errors relative to the intrinsic fluctuations in the chiral and continuum behaviour among the ensembles. Therefore we will focus on the TMR and the Hybrid method as implemented via Padé fits in the following.
Examples of our chiral and continuum extrapolations are shown in Fig. 3 while Table 7 contains an overview of results for the individual flavour contributions to a hvp µ at the physical point. We observe good agreement between the Hybrid and TMR methods. We also note that the inclusion of the finite-volume correction via the GS model produces a sizeable upward shift in (a hvp µ ) ud . This is also apparent from Fig. 4 Table 7: Summary of results for the hadronic vacuum polarization contribution (in units of 10 −10 ) at the physical point. The first error is statistical while the second denotes the systematic uncertainty as estimated via the variations listed in Table 6. The rightmost column contains the estimate for (a hvp µ ) ud including corrections for finite-size effects.
first concerns the impact of the uncertainty in the lattice scale: In order to make contact between the kernel function K(Q 2 ; m 2 µ ) and the VPFΠ(Q 2 ) computed on the lattice one must express the dimensionless momentum scale (aQ) in units of the muon mass. In our calculation the lattice spacing is known with a precision at the level of 1% (see Table 1). To assess the systematic error associated with scale setting we have repeated the chiral and continuum fits for the Hybrid method, using the upper and lower values of a as defined by the 1-σ bands. The variation of the lattice scale by ±1% increased the overall systematic error in (a hvp µ ) ud as estimated via the EFM by 1.8%. Given the ultimate precision goal of less than 1% uncertainty, this is a rather large systematic effect. For the TMR we have derived an entirely consistent estimate of the scale setting uncertainty using the representation of the kernel function K(x 0 ; m µ ). Details are presented in appendix B.2.
The second additional uncertainty is associated with the contributions from disconnected diagrams. In appendix D we present our calculation of quark-disconnected contributions on a subset of our ensembles (E5 and F6). The main result of that investigation is the derivation of a conservative upper bound on the magnitude of the disconnected contribution. Our findings indicate that quark-disconnected diagrams decrease the estimate of a hvp µ by at most 2%. As our final estimate for the hadronic vacuum polarization contribution we quote the result from the TMR including the finite-volume corrections based on the GS-parameterization. Adding the contributions from the light, strange and charm quarks we arrive at a hvp µ = (654 ± 32 stat ± 17 syst ± 10 scale ± 7 FV + 0 −10 disc ) · 10 −10 .
The quoted systematic error was estimated via the EFM considering the variations listed in the lower part of Table 6. The scale uncertainty (third error) amounts to the increase in the systematic error estimate when the lattice spacing is shifted by ±1 σ and the corresponding variations are included in the EFM procedure for the Hybrid method. As described in appendix C, we assign an uncertainty of 20% to the determination of the finite-volume shift in (a hvp µ ) ud . This produces an additional systematic error of ±7 · 10 −10 . Finally, we estimate that quark-disconnected diagrams reduce the value of a hvp µ by at most 10 · 10 −10 when the latter is computed using connected correlators only.
Our calculation has been performed in two-flavour QCD, and hence our results will be affected by the quenching of the strange and, to a lesser extent, the charm quark. Since we know of no reliable way of estimating the associated systematic effect, we leave it unspecified  The yellow vertical band denotes the result obtained from dispersion theory [3]. and caution the reader that this has to be taken into account when comparing our result to phenomenology or other lattice determinations. We add that our results are in good agreement with those of refs. [33,39] which were performed for N f = 2 + 1 + 1 flavours.

Conclusions
We have presented a lattice calculation of the hadronic vacuum polarization contribution to the muon g−2 addressing all sources of systematic error, except isospin breaking and the effects of dynamical strange and charm quarks. Given the overall uncertainty of 6% it is unlikely that our result, presented in eq. (37), is strongly biased by the omission of these effects. Our estimate is lower than the current value from dispersion theory but in agreement within the error of our calculation. Lattice determinations of a hvp µ have become more accurate in recent years, yet the target precision of < ∼ 1% has not been reached so far. While the statistical accuracy can be straightforwardly improved by an increased numerical effort, this is a lot more difficult for some of the various sources of systematics error.
In this paper we have investigated several complementary methods designed to control the infrared regime. One important lesson is the observation that this issue is strongly linked with the question of finite-volume effects. Our investigation of the long-distance regime of the vector correlator by means of the Gounaris-Sakurai parameterization of the pion form factor revealed that finite-volume effects are significant. They amount to a 5% shift in the value of a hvp µ for m π L ≈ 4 and near-physical pion masses. While this is consistent with similar estimates based on effective field theories (see, for instance, refs. [39,40,78,94]), a direct calculation, performed at sufficiently large m π L, which demonstrates that finite-volume effects are under control is still lacking. Based on the Gounaris-Sakurai model, we estimate that finite-volume effects are below the percent level when m π L > ∼ 6. Another important issue is the individual contribution from the charm quark, (a hvp µ ) c , which amounts to about 2% of the total value. Given that (a hvp µ ) c is quite sensitive to lattice artefacts, it is of vital importance to reliably control the continuum limit if one aims at sub-percent precision. Furthermore, scale setting has a large influence on the overall accuracy. Our analysis has shown that an extremely precise calibration of the lattice spacing -significantly below the percent level -is indispensable for a lattice determination of a hvp µ that is competitive with the dispersive approach.
Acknowledgments: The authors are indebted to Jeremy Green for the calculation of renormalization factors. We are grateful to our colleagues within the CLS initiative for sharing ensembles. Our calculations were partly performed on the HPC Clusters "Wilson" and "Clover" at the Institute for Nuclear Physics, University of Mainz. We thank Dalibor Djukanovic and Christian Seiwerth for technical support. We are grateful for computer time allocated to project HMZ21 on the BG/Q "JUQUEEN" computer at NIC, Jülich. This work was granted access to the HPC resources of the Gauss Center for Supercomputing at Forschungszentrum Jülich, Germany, made available within the

A Renormalization of the vector current
Here we describe the procedure used to determine the (mass-dependent) renormalization factor of the vector current from the quark-connected contribution to the three-point function and the two-point function where the operator O is given by O = ψ f γ 5 ψ f , and V loc µ,f is defined in eq. (17). Choosing the source-sink separation t s as t s = T /2 one can form the ratio as well as the difference d(t) ≡ R(t, T /2) − R(t + T /2, T /2).  By fitting d(t) to a constant Q V over a Euclidean time interval one can determine the renormalization factor Z computed on all ensemble used in this study. The renormalization condition of eq. (42) depends on the flavour f of the spectator quark. On ensemble E5 we have studied all possible combinations of f and f (i.e. ud, s and c). Our findings indicate that spectator quark effects are below 1%, with the strongest influence seen in the case of the renormalization of the charm quark contribution to the vector current.

B The QED kernel in the time-momentum representation
The vector correlator in the time-momentum representation is given in eq. (6). The master equation to compute a hvp µ from it is [51] a hvp with the momentum-space kernel given by 3

B.1 Derivation of a representation of the kernel function
Our goal is to obtain a simple and accurate representation off (t) which can be used straightforwardly in the expression for a hvp µ via eq. (43). Sincef (t) has units of GeV −2 and only involves the muon mass as an external scale, it is clear that m 2 µf (t) must be a dimensionless function in the variable (m µ t).
For the following derivation it is convenient to set the muon mass to unity and restore the units by dimensional analysis at the end of the calculation. The function f (ω 2 ) can be simplified (ω > 0), and hence f (ω 2 )/ω goes like 1/ω 2 at ω = 0.
The key observation is thatf (t) can be expressed in terms of the auxiliary functioñ Note that > 0 serves as an infrared regulator which is removed at the end of the calculation. In fact, we note that the regulation is only necessary for the first two terms in f (ω 2 ). One finds that the contribution of the second and third term in eq. (47) tog (t) can be expressed in terms of modified Bessel functions, K 0 and K 1 . The first term in eq. (47) is the most complicated: It involves the evaluation of the integral which satisfies The two linearly independent solutions of the homogeneous equation are e ± t . A particular solution I p (t) of the inhomogeneous equation can be found using the standard integral representation and the Laplace transform Realizing that can be set to zero, we arrive at the representation Noting that I p (0) = −1/4 and I p (0) = π/4, we impose the initial conditions and obtain the full solution up to terms of O( ), i.e.
The integral I p (t) can be expressed in terms of Meijer's G function [95]. In Mathematica [96], it can be evaluated by a built-in function Putting everything together, we havẽ From here one obtains straightforwardly, now restoring the units, where γ E = 0.57721566490153286061 . . . is Euler's constant. The expansion off (t) around the origin yields .
Note thatf (t) is not analytic at the origin, due to the appearance of terms proportional to ln(t) beyond fourth order. The expansion at large t yields m 2 µf (t) = 2π 2t2 − 4π 3t + 4π 2 (−1 + 4γ E + 4 ln(t)) + For a numerical evaluation, we propose the following. Up tot = 1.05, the expansion of eq. (58) around the origin provides an estimate off (t) with a relative accuracy better than 3.3 · 10 −6 . Beyond that point, the series can be used. Its accuracy is also better than 3.3·10 −6 for allt ≥ 1.05. Note that the integrand for a µ is expected to be very small beyond 4 fm, corresponding tot > 2.14; see Fig. 4 in [51].

B.2 Sensitivity of a hvp
µ to the lattice scale setting The representation for the kernel functionf derived above can be used to study the sensitivity of a hvp µ on the uncertainty in the determination of the lattice scale. Standard error propagation implies that the uncertainty ∆Λ on the observable Λ that sets the lattice scale translates into a corresponding uncertainty in a hvp µ according to where M µ ≡ m µ /Λ denotes the muon mass in units of Λ. To evaluate the derivative, we note that A short calculation then leads to As an example application, using the parameterization of the R-ratio in [51], which yields

C Finite-size effects in the time-momentum representation
In this appendix we address the finite-size effects on a hvp µ in the TMR and our ability to calculate them. Finite-size effects on the time-momentum correlator G ρρ (x 0 ) were computed in [54] based on the Lüscher formalism and the relation between the timelike pion form factor and finite-volume matrix elements [87,88]. Here we employ exactly the same method and therefore refer the reader to [54] for the relevant technical details. The goal of this appendix is to study the finite-size effects we expect on theoretical grounds at the simulation parameters used in the actual calculation presented in the main text. Several groups have studied finitesize effects on the hadronic vacuum polarization by theoretical means, see [78,94]. In any comparison, one must keep in mind that the finite-size effects depend on precisely which finite-volume representation of a hvp µ or the vacuum polarization one is using. We will compare our predictions quantitatively to the leading prediction of chiral perturbation theory.
The only input required in our analysis is the timelike pion form factor, including its phase, which coincides with the iso-vector p-wave ππ scattering phase. We use the phenomenologically successful Gounaris-Sakurai (GS, [92]) parameterization of the form factor  Table 9: Parameters of the Gounaris-Sakurai model used to explore finite-size effects on the various ensembles. P4 is a hypothetical ensemble at the physical pion mass. The width parameter at the physical pion mass is taken from [54], and is estimated from there for the other pion masses according as described in [54], noting that alternative parameterizations are available (see [97] and references therein). Clearly, the most important feature in the form factor is the ρ-resonance.
The main finite-size effect is that the finite-volume correlator falls off more rapidly than its infinite-volume counterpart, because the finite-volume spectrum is discrete and starts at a higher energy than 2m π . In order to proceed, we separate the correlator into two parts, t < t i and t > t i , with t i ≈ 1 fm. The reason for doing so is that the long-distance part can be analyzed using the lowlying energy-eigenstates on the torus. At shorter distances, the Poisson-resummed expression based on non-interacting pions should provide a good approximation to the finite-size effects for realistic m π L ≥ 4 [54]. As we show below, the finite-volume effects for the contribution to a µ from t < 1 fm are negligible for m π L ≥ 4 and m π < ∼ 300 MeV. Specifically, we define the short-and long-distance contributions computed on a finite torus as follows, Heref (t) is the QED kernel, given explicitly in appendix B. The Euclidean time t i represents the point beyond which the two-pion channel dominates the correlator. Using the Gounaris-Sakurai model combined with the Lüscher formalism for a > µ , as in [54], we obtain for the sets of parameters listed in Table 9 the estimates of the finite-size effects in Table 10. The effects are sizeable compared to the ultimate sub-percent accuracy goal. In addition to the lattice ensembles available to us, we also consider for illustration an ensemble at the physical pion mass and m π L = 4, labelled P4. For a < µ , we use the free-pion approximation to compute finite-size effects. Some details of this approximation are given in the next subsection.  in the TMR in units of 10 −10 , based on noninteracting pions for the 'short-distance' contribution a < µ and on the Gounaris-Sakurai model of the timelike pion form factor and the Lüscher formalism for the 'long-distance' contribution a > µ . The last column is discussed in section C.3. We used the values t i = (m π L/4) 2 /m π , t f = 1 fm and t cut = max(t i , 1.35 fm). The parameters used for the different ensembles are listed in Table 9.
We compute the finite-size effect from the part t < t i using eq. (67) and obtain the values quoted in Table 10, column 4. The small values indicate that the finite-size effects from the region below about 1 fm can be neglected for m π < ∼ 300 MeV and for m π L ≥ 4. If we compute the finite-size effect at large Euclidean times using free pions (using Eq. (66)), we obtain for instance 12.6 (P4, t i = 1.41 fm) 8.0 (G8, t i = 1.07 fm) We see that, although of the same order of magnitude as the finite-size effects in Table 10 (column 5) estimated using the Gounaris-Sakurai model in conjunction with the Lüscher formalism, the numbers in eq. (68) are smaller by a factor 1.5-2.0. For any fixed t, we expect the free-pion theory to predict the leading finite-size effect (O(e −mπL )) for L sufficiently large. However, at times t > 1 fm, many terms contribute significantly in the winding expansion eq. (67) at realistic parameters. It is then more expedient to use the sum over energy eigenstates as in eq. (66), however, with the energy levels and matrix elements taking into account ππ interactions via the Lüscher formalism. We conclude that the interactions between pions play an important role in estimating the finite-size effect in the t > 1 fm region at the typical volumes m π L ≈ 4. The Gounaris-Sakurai model also allows us to estimate a lower bound on the value of m π L for which finite-size effects in a hvp µ are below the level of 1%. From Table 10 we can read off that finite-size effects from the region t > 1.4 fm are as large as 3% for ensemble P4. By repeating the analysis for larger values of m π L we find that finite-size effects from the region t > 1.4 fm are reduced to about 1% when m π L ≈ 6. By contrast, finite-size effects from the region below 1.4 fm are already well below 1% for m π L = 4.

C.2 Reliability of the estimate of finite-size effects
To discuss the dependence of our theory estimate of the finite-size effect on the parameters, we focus on the ensemble G8, where the correction is sizeable. Using the GS model combined with the Lüscher formalism, we obtain at t = t i = 1.07 fm, while for free pions, we get for the same quantity 3.3 · 10 −10 . Thus at the turning point, where we switch from the free-pion to the interacting-pion case, the difference between the two predictions is moderate. This is a first indication that the overall prediction of the finite-size effect is not too sensitive to the turning point t i . Explicitly, we explore the dependence of the predicted finite-size effect on various parameters in Table 11. The result hardly changes under reasonable variations of t i , m ρ and Γ ρ . Of course the small observed variations do not reflect the full uncertainty due to the use of the Gounaris-Sakurai parameterization, the corrections to the finite-size effect at t < t i due to pion interactions and internal structure, etc. We think that the genuine finite-size effects on a hvp µ (i.e. the sum of column 4 and 5 in Table 10) are correctly estimated to within 20% in our approach.
We have also performed a sanity check by comparing our prediction for finite-size effects to the direct lattice QCD data in [39], where at one set of quark masses, results for a hvp µ at three volumes are available: within the uncertainties, our estimate for the volume-dependence of dΠ dQ 2 | Q 2 =0 is fully consistent with the numerical data. In the comparison, we assume that finite-size effects are dominated by the iso-vector contribution to a hvp µ , since the iso-scalar ω and φ resonances are extremely narrow.

C.3 Single-exponential extension of the time-momentum correlator
Since in practice an extension of the vector correlator is used at long distances, we introduce where t cut > t i is the point beyond which the one-exponential extrapolation of the finitevolume correlator G8: parameter varied 10 10 (a µ (∞) − a µ (L)) t i = 1.   Table 9; they lead to a µ (∞) − a µ (L) = 16.0 · 10 −10 (sum of column 4 and 5 in Table 10).
is used, based on the effective mass and amplitude determined at time t f ; explicitly, The reason for considering a >,xpol µ (t i , t f , t cut , L) is that due to the deteriorating signal-to-noise ratio on the vector correlator at large distances, some form of extrapolation is required in practice to be able to integrate to t = ∞.
We indicate in the last column of Table 10 what error one incurs by replacing the correlator by its one-exponential extension beyond t cut . As compared to the genuine finite-size effect (column 5 of the table), the additional systematic error is relatively modest until one reaches the ensembles with m π below 200 MeV. At this point, the result is also quite sensitive to the time t f where the effective mass is determined. On ensemble G8 for instance, we obtain Thus for ensembles with m π < ∼ 200 MeV, the single-exponential extension is clearly inadequate once the precision goal on a hvp µ is 5% or better.

C.4 Uncertainty in the determination of the ρ-mass and decay width
In the absence of a full dedicated study of the spectroscopy in the iso-vector vector channel, in section 4.2 we have assumed the GS form of the timelike pion form factor and used a simplified procedure to determine the parameters (m ρ , Γ ρ ) of the model. On our ensemble G8 with the lightest pion mass, we assumed that the ground state had an energy of E 0 = 2 m 2 π + (2π/L) 2 corresponding to non-interacting pions in a p-wave, while the energy of the first excited state was identified with the parameter m ρ of the GS model. We have investigated how reliable these assumptions are using the GS model; see Fig. 5. Especially the first excited state corresponds to the ρ-mass to sub-percent accuracy for a wide range of parameters. The deviation of the ground state from the non-interacting-pions predictions is at the 3-4% level. At our present level of accuracy, this is a sufficient level of control to avoid a significant bias in the determination of the first excited state, since the ground state contributes with a relatively weak amplitude to the vector correlator.  Figure 5: Corrections to energy levels relative to the naive expectation of a non-interacting, p-wave two-pion state and a ρ-state, for parameters corresponding to ensemble G8 and assuming the GS pion form factor. Left: correction to the expectation E 0 = 2 m 2 π + (2π/L) 2 for the ground-state energy as a function of the width Γ ρ , for three values of the mass m ρ . Right: correction to the expectation E 1 = m ρ .

D Determination of the quark-disconnected contribution
In this appendix we provide the details of our calculation of the quark-disconnected contribution to a hvp µ , which has been performed using the TMR formulation (see also contribution 2.16 in [98]). Analytic analyses of disconnected contributions have been presented in [99,100]. For our discussion it is useful to recall the expression for a hvp µ in the TMR, i.e.
where K(x 0 ; m µ ) is defined in eq. (8). In the following we restrict the analysis to the contributions from the u, d and s quarks only, so that the electromagnetic current is given by After performing the Wick contractions one can identify the connected and disconnected parts as where G ud and G s are defined according to eq. (22), and the total disconnected contribution G disc (x 0 ) is given by The superscripts indicate whether the contribution involves only light (ud), strange (s) or both (ud, s) quark flavours (note that we work in the isospin limit, m u = m d ).
In ref. [35] it was shown that G disc (x 0 ) factorizes according to  Table 12: Details of the evaluation of quark-disconnected contribution G disc (x 0 ) (see eq. (76)). N r denotes the number of stochastic sources per timeslice, while x * 0 represents the Euclidean time at which the ratio G disc (x 0 )/C ρρ (x 0 ) is replaced by its asymptotic value. The upper bound on the size of the quark-disconnected contribution to a hvp µ is given by ∆a hvp µ .
where ∆ f (x 0 ) for f = (ud), s is given by and S f denotes the quark propagator of flavour f . Statistically accurate results for quantities such as ∆ f require "all-to-all" propagators which are commonly computed using stochastic noise sources. In [35] it was shown that the statistical accuracy of G disc (x 0 ) can be significantly enhanced when ∆ ud and ∆ s are computed using the same random noise vectors, since the correlations between the light and strange quark contributions largely cancel the stochastic noise.
In our determination of G disc (x 0 ) we have used stochastic sources in conjunction with a hopping parameter expansion (HPE) of the quark propagator [101], suitably adapted to the case of O(a) improved Wilson quarks [102]. The calculation was performed at our intermediate value of the lattice spacing at pion masses of 437 and 311 MeV, respectively (ensembles E5 and F6). The all-to-all propagators for the light and strange quarks were computed by employing a 6th order HPE in combination with N r stochastic U(1) noise vectors η k ( x), k = 1, . . . , N r on each timeslice. Further details are listed in Table 12.
Results for G disc (x 0 ) on the two ensembles under study are shown in Fig. 6. While a small but non-zero signal is observed for x 0 /a < ∼ 8 the disconnected contribution G disc (x 0 ) vanishes within errors for larger values of x 0 . At small times the disconnected contribution is only about 0.005% of the connected one, and hence we conclude that the vector correlator G(x 0 ) is completely dominated by the connected part in the region x 0 < ∼ 0.5 fm. The fact that the disconnected contribution is small where it can be resolved does not, however, imply that it is negligible. Using our data we can derive an upper bound on the error which arises if one were to neglect the disconnected contribution altogether. To this end it is useful to recall the isospin decomposition of the electromagnetic current shown in eq. (13), which gives rise to the iso-vector (I = 1) correlator G ρρ and its iso-scalar counterpart G I=0 (see eq. (15)). The iso-vector correlator G ρρ (x 0 ) contains only quark-connected diagrams; it is related to the connected light quark contribution G ud (x 0 ) via G ρρ (x 0 ) = 9 10 G ud (x 0 ). of ensemble F6. Thus, there is no visible trend for distances x 0 < ∼ 1.7 fm that the ratio approaches its asymptotic value of −1/9. In order to derive a conservative upper bound on the quark-disconnected contribution we assume that the ratio of eq. (82) drops to −1/9 at the time x * 0 where the accuracy of the data is insufficient to distinguish between zero and the expected asymptotic value. In other words, we set If we write the hadronic vacuum polarization contribution a hvp µ as the sum of the quarkconnected and -disconnected contributions, a hvp µ = (a hvp µ ) con + (a hvp µ ) disc , we can define ∆a hvp µ := (a hvp µ ) con − a hvp µ (a hvp µ ) con ≡ − (a hvp µ ) disc (a hvp µ ) con , which is the relative size of the disconnected and connected contributions, and (a hvp µ ) disc is given by After inserting eqs. (86) and (80) we obtain the maximum estimate of the quark-disconnected contribution as The resulting estimates for the relative contribution ∆a hvp µ are listed in Table 12.