The leading hadronic contribution to $(g-2)_\mu$ from lattice QCD with $N_{\rm f}=2+1$ flavours of O($a$) improved Wilson quarks

The comparison of the theoretical and experimental determinations of the anomalous magnetic moment of the muon $(g-2)_\mu$ constitutes one of the strongest tests of the Standard Model at low energies. In this article, we compute the leading hadronic contribution to $(g-2)_\mu$ using lattice QCD simulations employing Wilson quarks. Gauge field ensembles at four different lattice spacings and several values of the pion mass down to its physical value are used. We apply the O($a$) improvement programme with two discretizations of the vector current to better constrain the approach to the continuum limit. The electromagnetic current correlators are computed in the time-momentum representation. In addition, we perform auxiliary calculations of the pion form factor at timelike momenta in order to better constrain the tail of the isovector correlator and to correct its dominant finite-size effect. For the numerically dominant light-quark contribution, we have rescaled the lepton mass by the pion decay constant computed on each lattice ensemble. We perform a combined chiral and continuum extrapolation to the physical point, and our final result is $ a_\mu^{\rm hvp}=(720.0\pm12.4_{\rm stat}\,\pm9.9_{\rm syst})\cdot10^{-10}$. It contains the contributions of quark-disconnected diagrams, and the systematic error has been enlarged to account for the missing isospin-breaking effects.


I. INTRODUCTION
Electrons and muons carry a magnetic moment, which is correctly predicted by Dirac's original theory of the electron to within a permille of precision. The proportionality factor between the spin and the magnetic moment of the lepton is parameterized by the gyromagnetic ratio g. In Dirac's theory, g = 2, and one characterizes the deviation of g from this reference value by a = (g − 2) /2. Testing the ability of Quantum Electrodynamics (QED) to correctly predict this precision observable has played a crucial role in the development of quantum field theory in general. Presently, the achieved experimental precision of 540 ppb on the measurement of the anomalous magnetic moment of the muon [1], a µ , requires the effects of all three interactions of the Standard Model (SM) of particle physics to be included in the theory prediction. In fact, a tension of about 3.5 standard deviations exists between the SM prediction and the experimental measurement. For reviews on the subject, we refer the reader to [2][3][4].
Presently, the E989 experiment at Fermilab is performing a new direct measurement of a µ [5], and a further experiment using a different experimental technique is planned at J-PARC [6]. The final goal of these experiments is to reduce the uncertainty on a µ by a factor of four. A reduction of the theory error is thus of paramount importance, as the first results from the Fermilab experiment are expected within the next few months. These will likely reach the same precision as the current world average.
On the theory side, the precision of the SM prediction for a µ is completely dominated by hadronic uncertainties. The leading hadronic contribution enters at second order in the finestructure constant α via the vacuum polarization and must be determined at the few-permille level in order to match the upcoming precision of the direct measurements of a µ . In this paper we undertake a first-principles lattice QCD calculation of this hadronic contribution (see [7] for a recent review of previous lattice results). A further hadronic effect, the lightby-light scattering contribution which enters at third order in the fine-structure constant, currently contributes at a comparable level to the theory uncertainty budget and is being addressed both by dispersive and lattice methods (see [8][9][10] and references therein).
Our calculation of the hadronic vacuum polarization to the anomalous magnetic moment of the muon, a hvp µ , fully includes the effects of the up, down and strange quarks, while the charm quark (whose contribution to a hvp µ is small) is treated only at the valence level. We use ensembles of SU(3) gauge field configurations generated with an O(a) improved Wilson quark action as part of the Coordinated Lattice Simulations (CLS) initiative [11,12]. In particular, the generation of a physical-mass ensemble [13] (labelled E250) was largely motivated by the goal of improving the lattice determination of a hvp µ . We use four different lattice spacings to control the continuum limit, and the (u, d, s) quark masses are varied at constant average quark mass in order to perform a chiral interpolation to the physical values of the quark masses [11]. Our calculation is performed at equal up and down quark masses, and no QED effects are included; however, in the future both of these isospin-breaking effects will be taken into account as corrections [14,15].
Lattice QCD, which is formulated in Euclidean space, is well suited for computing a hvp µ , since the latter only involves the two-point function of the hadronic component of the electromagnetic current at spacelike momenta [16]. In this work we employ the representation of a hvp µ as a Euclidean-time integral over the two-point function in the time-momentum representation (TMR) [17], i.e. projected to vanishing spatial momentum. This representation does not require a parameterization of the vacuum polarization function and has a clear spectral interpretation in terms of vector hadronic states in the center-of-mass frame. The main difficulty in obtaining a hvp µ with good statistical precision is that it probes the TMR correlator at Euclidean times well beyond 2 fm, where its relative precision deteriorates rapidly. Therefore a dedicated treatment of the tail of the correlator which does not compromise the first-principles nature of the calculation is needed. Here the spectral representation of the correlator plays a central role.
An important source of systematic uncertainty is the correction to a hvp µ due to the use of a finite spatial torus. On our lattice ensembles, this finite-size effect (FSE) mostly stems from the tail of the isovector component of the TMR correlator. Thanks to precise relations [18,19] between the properties of the discrete quantum states on the torus and the pion form factor at timelike momenta, we are able to correct for the dominant part of the FSE. Finally, the quark-disconnected diagrams, while making only a few-percent contribution to a hvp µ , require a dedicated set of calculations for their evaluation, which demand a large computingtime investment.
The rest of this paper is organized as follows. Section II describes the methodology followed in our calculation, including the renormalization and improvement of the TMR correlator and the treatment of the charm contribution. Section III presents our lattice data and the extraction of the observable a hvp µ on each individual lattice ensemble. In section IV, the lattice-spacing and quark-mass dependence of these intermediate results is fitted in order to arrive at our final result. Finally, we compare the latter with phenomenological as well as other recent lattice determinations in section V.

A. Time-momentum correlators
We start by providing all relevant relations in the continuum and infinite-volume Euclidean theory. In the time-momentum representation (TMR), the leading-order hadronic vacuum polarization contribution to (g − 2) µ is given by the convolution integral where an analytic expression for the QED kernel function K(t) is given in Appendix B of Ref. [20], and is the spatially summed QCD two-point function of the electromagnetic current J = 2 3ū γu− 1 3d γd − 1 3s γs + 2 3c γc. In isospin-symmetric QCD, we can write where G f (t) denotes a quark-connected contribution associated with flavour f and G disc (t) is the quark-disconnected contribution. An alternative decomposition based on the isospin quantum number I yields Physically, the latter decomposition is more transparent. In particular, at light pion masses the dominant finite-size effects, as well as a logarithmic singularity as m π → 0, only concern the isovector contribution, a hvp,I=1 µ . Computationally however, the disconnected contributions are obtained very differently from the connected ones: they are costly and amount only to a few percent of the total. Therefore, in our numerical analysis there is an interesting interplay between the two choices of bases to compute a hvp µ . With m µ the muon mass, the kernel behaves as K(t) ∼ π 2 9 m 2 µ t 4 for t m −1 µ and as K(t) ∼ 2π 2 t 2 for t m −1 µ . Since the lattice data for the correlator G(t) is in lattice units, the muon mass must be known in those units, am µ . The knowledge of the lattice spacing in GeV −1 thus plays a crucial role in a precision determination of a hvp µ [20,21]. There are then two ways to proceed. In lattice QCD, where often the physical quark masses are reached only after an extrapolation or interpolation, a hvp µ can either be calculated using the fixed, physical value of m µ = 105.66 MeV; or the muon mass can be rescaled by a quantity with dimension of mass known experimentally [22]. In our calculation, we have explored both paths. In our final results, we adopt the "rescaling strategy" for the connected light contribution. As a rescaling quantity, we choose the pion decay constant f π , so that we set 1 on every lattice ensemble. Our choice is motivated, first, by f π being determined precisely and reliably, both in phenomenology and on the lattice; and secondly, since f π increases with the pion mass, this choice has the effect of making the m π dependence of a hvp µ weaker. To intuitively understand the effect of the rescaling, it is instructive to consider the calculation of the anomalous magnetic moment of the electron; in this case, obtaining a hvp e = (4α 2 /3)m 2 e Π 1 requires computing the time moment Π 1 ≡ (1/12) ∞ 0 dt t 4 G(t). Thus the rescaling simply amounts to computing the dimensionless quantity f 2 π Π 1 , and converting the result into a hvp e by using the phenomenological value of (m 2 e /f 2 π ).

B. Simulation parameters
Our work is based on a subset of the Coordinated Lattice Simulations (CLS) ensembles with N f = 2 + 1 dynamical quarks. They are generated [11] using the open-QCD suite 2 [23] and are based on the O(a)-improved Wilson-Clover action for fermions, with the parameter c sw determined non perturbatively in Ref. [24], and the tree-level O(a 2 ) improved Lüscher-Weisz gauge action. The ensembles used in this analysis were generated at a constant value of the average bare quark mass such that the improved bare couplingg 0 is kept constant along the chiral trajectory [11]. In particular, five of the ensembles are at the SU(3)-symmetric point, m u = m d = m s . The parameters of the simulations are summarized in Table I.
Results are obtained at four values of the lattice spacing in the range a = 0.050−0.086 fm. The scale setting was performed in Ref. [25] using a linear combination of the pion and kaon decay constants with a precision of 1%. The pion masses used in our determination of a hvp µ lie in the range m π ≈ 130 − 420 MeV. All the ensembles included in the final analysis satisfy m π L > 4. Furthermore, at two values of the pion mass (m π = 280 and 420 MeV), two ensembles with the same bare lattice parameters but different volumes are used to study finite-size effects. These ensembles with smaller volumes are not included in the final analysis and are marked by an asterisk in Table I.
All ensembles have periodic boundary conditions (BC) in space. In the time direction, ensembles E250 and B450 have periodic BCs, while all others have open temporal BCs. The choice of open boundary conditions was made in order to address the issue of long autocorrelation times associated with the topological charge at small lattice spacing [26]. Our use of ensembles with open BCs constitutes part of our motivation for employing correlators in the time-momentum representation. The boundary couples to a tower of states with vacuum quantum numbers. Therefore, in order to extract vacuum correlators, sources and sinks of correlation functions should be placed at a sufficient Euclidean-time separation away from the boundaries 3 . On the ensembles with periodic temporal BCs on the other hand, we exploit the translation invariance in time to increase statistics.
For all ensembles, except E250, the TMR correlation functions are computed using point sources, randomly distributed in space and in the center of the lattice in the time direction. As described in the next subsection, we use the local vector current at the source and both the local and the conserved vector currents at the sink. For the ensemble E250, propagators are estimated using stochastic sources, with noise partitioning in spin, colour and time [27,28]. Each source has support on a single, randomly chosen timeslice. To improve statistics, the TMR correlator in Eq. (2) is averaged over the three spatial directions. Errors are estimated throughout the calculation using the jackknife procedure with blocking in order to take into account auto-correlation effects.
In addition to the direct calculation of the TMR correlators, the auxiliary calculation of the ππ I = = 1 scattering phase plays an important role in our determination of a hvp µ . In Ref. [29], it has been determined on ensembles C101, N401, N200, D200 and J303. On all these ensembles except C101, the pion form factor at timelike kinematics has also been determined in [29]. As compared to the latter reference, the number of gauge configurations used for our spectroscopy calculation on ensemble D200 has roughly been doubled. Additionally, we have performed a spectroscopy calculation on ensemble N203 with a statistics of about 200 gauge configurations.
We have computed the quark-disconnected contribution to a hvp µ on ensembles N401, N203, N200, D200 and N302. This selection provides us with a handle on the discretization effects at m π 345 MeV and m π 285 MeV, and allows us to investigate the chiral behaviour of the disconnected contribution via the fixed lattice-spacing sequence of ensembles N203, N200, D200. The disconnected quark loops are computed using four-dimensional, hierarchically probed noise sources [30] with 512 Hadamard vectors. More technical details on our implementation can be found in [31].

C. Lattice correlators, renormalization and O(a) improvement
There are two commonly used discretizations of the vector current in Wilson lattice QCD, the local and the conserved current. For a single quark flavour q, their expressions are In our calculation of correlation functions, we always place the local vector current at the origin in Eq. (2); at point x, we use either the local or the conserved vector current. This provides us with two discretizations of the TMR correlator which share the same continuum limit. The conserved vector current has the advantage of not undergoing any renormalization or flavour-mixing. As for the flavour structure, we note that the electromagnetic current can be decomposed in the SU(3) Gell-Mann basis as where V a µ =ψγ µ λ a 2 ψ, withψ = (ū,d,s). Therefore, the local current only requires the non-singlet renormalization factor Z V . The charm-quark contribution is treated separately, at the "partially quenched" level; our treatment of this (small) contribution is described in the next subsection.
We have implemented the Symanzik O(a) improvement programme as described in Ref. [32]. Since our lattice action is O(a) improved, we now describe the improvement and renormalization of the vector currents in order to consistently carry out the Symanzik programme. The first step is to add to the local vector current an additive O(a) counterterm with a tuned coefficient c L V (respectively c C V for the conserved current) compensating chiral-symmetry violating effects in on-shell correlation functions, where ∂ ν denotes the symmetric lattice derivative 4 and Σ a µν (x) = − 1 2ψ (x)[γ µ , γ ν ] λ a 2 ψ(x). The second step, which is only required for the local current, is to take into account the following renormalization pattern, We denote by V L,0 µ = 1 2ψ γ µ ψ the flavour-singlet current, and the mass-dependent renormalization factors are given by [33,34] Here (m q,l , m q,l , m q,s ) are the bare subtracted quark masses, m av q is their average andg 0 is the O(a) improved bare coupling. We note that the mixing coefficient Z 80 is of order a and vanishes in the SU(3)-flavour-symmetric limit. We use the values of the renormalization factor Z V , the critical hopping parameter κ crit , as well as the improvement coefficients b V , b V , c L V and c C V , which are functions of the bare coupling g 0 , determined recently in [34]. There it was shown that the obtained values of Z V differ by percent-level O(a 2 ) effects from an independent high-precision determination [35]. The improvement coefficient f V (g 0 ), which is of order g 6 0 and only affects the isoscalar contribution to a hvp µ , is neglected. We estimate that the systematic error incurred by this approximation is at present negligible. Strictly speaking, the connected strange correlator taken in isolation requires an independent, partially quenched improvement coefficient in the mass-dependent part of the renormalization factor in order to be consistent with O(a) improvement; however, we have neglected this effect.
The desired quantity a hvp µ is obtained using Eq. (1), where the integral is replaced by a sum over timeslices. Note that in the improvement terms entering the TMR correlator, only the temporal derivative of the tensor current contributes. For the connected light and strange contributions, we have compared the use of the symmetric lattice derivative in Eq. (8) with an alternative implementation where an integration by parts is used in order to apply the temporal derivative on the QED kernel K(t), and found the difference to be negligible; therefore we have used the symmetric lattice derivative throughout. For the charm contribution, however, we have found it advantageous to use the discrete derivative on the 'away' side, i.e. in such a way that the vector-tensor correlator is not evaluated at a shorter time separation than the vector-vector correlator itself [36]. Finally, we remark that we do not include the O(a 2 ) term consisting of the correlation of two tensor currents.

D. Treatment of the charm contribution
We treat the charm quark at the partially quenched level: it does not appear in the simulated action, nor do we include the contribution of quark-disconnected diagrams containing charm loops. Given that the charm contribution is about two percent of the total, these approximations appear fully sufficient at our present level of precision.
The first task is to tune the value of the bare charm quark mass on each lattice ensemble. The mass of the ground state pseudoscalar cs meson is computed for several values of κ c , using stochastic sources with colour, spin and time dilution. The value of κ c used in the calculation of a hvp µ is then obtained from a linear interpolation of the squared mass of the lightest cs meson in 1/κ c to the point where this mass equals the experimental value of the D s meson mass.
We perform a dedicated determination of the multiplicative, mass-dependent renormalization factor Z c V for the local charm current on every lattice ensemble. The determination is based on requiring the charm quantum number of the pseudoscalar cs meson to be exactly unity. It follows the method used in [34] for the light isovector current, and a similar method was already used in [20]. As for the improvement coefficients c L V and c C V , we use the same values as for the u, d, s quark flavours. The results for κ c and Z c V are given in Table IV, while the individual pseudoscalar cs meson masses used for the determination of κ c are collected for reference in Table IX of Appendix B.
E. Infrared aspects of a hvp µ : correlator tails, finite-size effects and the chiral limit There are a number of aspects of the calculation of a hvp µ related to the long-distance physics of vector correlators that are best discussed together. Here, we summarize our understanding of these issues before applying it to the treatment of lattice data.
In preparation, recall that the TMR correlator can be written, via the spectral decomposition in finite volume, as the sum of the (positive) contributions of individual vector states. In particular, only isovector vector states contribute in the correlator where the amplitudes Z n are real, and the discrete, ordered energies E n are real and positive. A similar expression holds for the isoscalar correlator G I=0 (t).

Controlling the long-time tail of the TMR correlators
The contribution of the tail of the correlator to a hvp µ is enhanced by the QED kernel. Yet the correlator is affected by a growing statistical error, as well as a large relative finite-size effect. We discuss these two issues in turn.
In order to handle the tail of the correlators, two types of treatment have been proposed. Both are based on the fact that at large Euclidean times, a few terms in the sum of Eq. (13) saturate the correlator to a high degree of precision, which was one of the motivations for introducing the time-momentum representation [17]. In the first type of treatment, one explicitly constructs an extension of the correlator for t > t c , motivated by the spectral representation (13). The simplest incarnation of this method, partly used in our earlier calculation [20], is to keep only the lightest of those states and thus to perform a one-exponential fit to the correlator for Euclidean times around t c . When a dedicated spectroscopy calculation is available, several energy levels E n as well as the overlaps Z n can be used, so that the summed contributions of these states already saturate the TMR correlator at smaller Euclidean times.
A second type of treatment consists in bounding the Euclidean correlator from above and below [37][38][39], exploiting the positivity of the prefactors Z 2 n /(2E n ), where N = 0 in the simplest variant, and E eff (t) ≡ − d dt log G(t) is the "effective mass" of the correlator. As a refined variant of this method, a dedicated spectroscopy calculation delivering the energies and matrix elements of the N lowest-lying states allows one to improve the control over the tail by applying the bound Eq. (14) to the subtracted correlator A challenge one eventually faces in exploiting lattice spectroscopy information is that the number of states required to saturate the TMR correlator at a given t c increases with decreasing pion masses and (roughly proportionally) with the volume. However, for the ensembles used in this work, the number of states needed is at most four.
Finite-size effects on a hvp µ in the time-momentum representation We now come to the closely related issue of the finite-size effect on the observable a hvp µ calculated in the time-momentum representation. At asymptotically large volumes, the finite-size effect is of order e −mπL and can be computed in chiral perturbation theory [40][41][42]. At low pion masses, the leading finite-size effect is expected to come from the ππ channel, and thus affects the isovector channel only, G I=1 (t). Working in the flavour decomposition of Eq. (3), we take this observation into account by applying 10/9 of the isovector finite-size correction to the connected light-quark contribution, and −1/9 of the same correction to the disconnected contribution.
Looking at the finite-size effect on the correlator as a function of Euclidean time, it has been pointed out [17,41] that for a given spatial box size L, the tail of the correlator is affected by an unsuppressed finite-size effect. One may define a time t i beyond which the finite-size effect becomes sizeable. While t i grows with L, we find that the overall finite-size effect on a hvp µ is dominated by the tail in our present calculation. For m π L = 4 − 5, the tail of the finite-volume isovector correlator is accurately described by the contribution of a handful of energy eigenstates; this point will be illustrated in Fig.  2. On the other hand, the tail of the infinite-volume correlator can be obtained from the timelike pion factor. Thus, knowledge of this form factor allows one to correct the tail of the isovector correlator [17]. In this work, we apply the same finite-size correction method as in our previous calculation [20], parameterizing the pion form factor with the Gounaris-Sakurai (GS) model [43]. While too simplistic a model for a study of the form factor for its own sake [29,44], we expect it to be sufficient for the purpose of reducing the residual finite-size effects to a level that is small compared to our current statistical precision. We emphasize that we only use the GS parametrization of the pion form factor for the finite-size correction, and not for the treatment of the tail of the correlators.
The chiral dependence of a hvp µ The TMR correlator for non-interacting pions was given in Ref. [41]. For massless pions, it is given by G(t) = 1/(24π 2 |t| 3 ); combined with the asymptotic form of the QED kernel for a finite muon mass, K(t) ∼ 2π 2 t 2 , this contribution generates a logarithmic divergence, which is made finite by a small but finite pion mass and then yields This result and further terms in the expansion have been derived in [45], where the systematics of the chiral extrapolation has been studied in detail. The asymptotic form (16) only becomes a decent approximation for m π /m µ well below 1/10. Thus this logarithmic divergence is largely irrelevant when describing the pion-mass dependence of a hvp µ in the range 130 < m π /MeV < 300. On the other hand, if m µ m π and both are small compared to the ρ meson mass, one finds the leading behaviour It turns out that this asymptotic form is rather robust, holding down to fairly small values of m π /m µ . In fact, within the framework of chiral perturbation theory at next-to-leading order 5 underlying Eqs. (16) and (17) At physical quark masses, the overall magnitude of expression (17) is enhanced by the (squared) pion form factor at timelike kinematics. In addition, the contribution of the ππ states with a center-of-mass energy well below the ρ-meson mass is numerically subdominant compared to the resonant contribution. The ρ-meson mass depends only mildly on the light-quark mass, and thus the steep behaviour predicted by Eq. (17) as a function of m π is superimposed on a larger, more slowly varying contribution. In our chiral extrapolations, presented in section IV, we use these observations to construct suitable fit ansätze for the chiral extrapolation.
The singular chiral behaviour comes from the isovector channel, while we expect the isoscalar channel to have a much milder dependence on the pion mass. Working in the basis of Eq. (3), the singular chiral behaviour is split between the connected light-quark contribution and the disconnected contribution. Indeed, in the limit that m µ and m π are much smaller than the hadronic scale, we have a hvp, disc µ = − 1 9 a hvp µ , and hence, from Eq. (17), For orientation, we note that if one inserts the physical pion mass into this expression, one obtains a hvp, disc µ = −10 × 10 −10 , and we expect this value to be further enhanced by the pion form factor. The important point is that the singular chiral behaviour present in the connected light-quark contribution to a hvp µ must be present in the disconnected contribution as well, with a relative factor of −1/10.

III. RESULTS
In this section we describe the main features of the TMR correlators obtained on the different lattice ensembles with a view to computing a hvp µ . Particular attention is devoted to the correlators at Euclidean times in the range [1.5, 4.0] fm. In the rescaling of the muon mass, we use the values of af π values given in Table VI, corrected for finite-size effects [48] and interpolated via a global fit in the pion mass and the lattice spacing.

A. The quark-connected contributions
The integrand of Eq. (1) for the connected light, strange and charm contributions is displayed in Fig. 1 for our two ensembles with quark masses closest to their physical values. The left (right) panel corresponds to a pion mass of about 200 MeV (131 MeV). The light contribution is clearly very dominant; note that the charm and strange contributions have been scaled by a factor of six for better visibility. On a given ensemble, the integrand peaks at increasingly longer distances as one goes from the charm to the strange to the light quarks, and the tail becomes more extended. At the same time, the statistical precision deteriorates. Comparing the left to the right panel, it is clear that the light contribution becomes harder to determine with the desired precision as the physical quark masses are approached. Nevertheless, these plots by themselves do not fully reflect all the known constraints on the TMR correlator, which is well known to be given by a sum of decaying exponentials with positive coefficients, as discussed in section II E.
Having described the state-of-the-art methods to handle the tail of the correlation function in section II E, we now describe how we applied these methods to our data. For the strange and charm quark contributions, the TMR correlator is determined so accurately that practically no particular treatment of the tail is needed. We apply the bounding method, Eq. 14 with N = 0, and obtain the results given in Table IV.
As for the connected contribution of the light quarks, our choice for the final analysis is again the bounding method on all ensembles; the only exception is the physical-pion-mass ensemble E250, to which we return below. In applying Eq. (14), we employ the expression containing the effective mass as a lower bound, and use as an estimate for the lowest-lying energy level in the channel the energy obtained by a one-exponential fit to the tail of the TMR correlator. On ensemble D200, on which the ground state lies clearly below the ρ mass and has a relatively weak coupling to the vector current, we use the auxiliary spectroscopy calculation to determine its energy. We find it to be close to, but slightly below the value  corresponding to two non-interacting pions, E free Table VI contains our results for the connected contributions of the light quarks.
As discussed in detail in the next subsection, the improved statistical precision gained by exploiting spectroscopic information can be quite significant for light pion masses, m π 200 MeV. Indeed we find that on the physical mass ensemble E250, on which we do not have direct spectroscopic information, we cannot achieve a comparable control over the statistical and systematic error with the simplest variant of the bounding method. Therefore we proceed as follows. The isovector vector energy levels computed on ensembles N203, N200 and D200 allow us to determine the scattering phase in the I = = 1 ππ channel [29] for energies up to the four-pion threshold via the Lüscher formalism [18] 6 . The scattering phase is well described by the effective range formula, with k ≡ 1 2 E 2 − 4m 2 π and k ρ being the value of k for E = m ρ . The parameters m ρ and Γ ρ correspond to the ρ meson mass and width. Furthermore, it has been observed in lattice simulations that parameterizing the width by the coupling g ρππ only has a weak pion-mass dependence. Therefore, we extrapolate the parameters (m ρ , g ρππ ) determined on the ensembles N203, N200 and D200 (see Table VIII) to obtain their values for the pion mass corresponding to ensemble E250. Using these values, we can predict the low-lying energy levels E n on ensemble E250 by using the Lüscher correspondence between them and the scattering phase in reverse. In order to obtain an extension of the TMR correlator on E250, we then fit the squared amplitudes Z 2 n , given the energy levels. Note that this can be formulated as a linear fit.
In our final choice of parameters, we fit the TMR correlator on E250 in the interval 26 < t/a < 37. Then the TMR is summed from t = 0 to t = 28a and the multi-exponential extension is used beyond that time. The numbers given for E250 in Table VI are  On ensemble D200 at m π = 200 MeV, we have detailed information on the scattering phase and the timelike pion form factor. We can thus test the validity of the procedure we applied on the physical pion-mass ensemble E250, described in the previous subsection.
Thus on D200 we apply and compare four different methods to handle the tail of the light connected correlator: (1) the bounding method without subtractions (N = 0); (2) the bounding method after subtracting the contribution of N = 2 states; (3) the extension of the correlator using the auxiliary information on the first two energy levels E n and their amplitudes Z n ; (4) the extension of the correlator using the auxiliary information on the first two energy levels E n , but fitting the amplitudes to the TMR correlator.
One motivation for comparing these particular methods is that on E250, we cannot apply the second or third method, while the first method would result in a large statistical error. Therefore, we apply the last method on E250, and presently test whether it gives consistent results on ensemble D200. µ with the f π -rescaled muon mass using the extension of the connected light (local-local) correlator using N = 2 energy levels on ensemble D200. On the left (method 3), the amplitudes corresponding to energy levels were predetermined in a spectroscopy calculation, while on the right (method 4), they are fitted to the TMR correlator. Results based on 1100 gauge configurations.   (1) and (2), as a function of the time t c at which the upper and lower bounds start to be used instead of the TMR correlator itself. The values are consistent with each other, however method (2) yields a significantly reduced statistical uncertainty. This outcome is not surprising, since important auxiliary information is used in method (2).
A comparison of methods (3) and (4) is shown in Fig. 4, showing the resulting a hvp µ as a function of the time t c at which the TMR correlator is replaced by the multi-exponential extension. The result of method (4) is consistent with that of method (3), albeit with an enlarged statistical uncertainty. In addition we have checked that the values of the amplitudes of the first two states as extracted from the fit in method (4) are well consistent with their direct spectroscopic determination. Table II presents the results obtained on D200 with the four different methods.

C. Finite-volume effects
As explained in section II E, in the isospin basis, we would correct the I = 1 correlator for finite-size effects stemming from the ππ states, and neglect such effects on the I = 0 correlator. However, we work in the basis of Eq. (3). In this basis, such a correction  corresponds to applying an additive finite-size correction to the connected light contribution ( 5 9 G l (t)), weighted by a factor of 10/9 relative to the correction of the I = 1 correlator. At the same time, the disconnected contribution G disc must be corrected by −1/9 of the I = 1 correction. It is indeed well known that the tail of G disc (t) is given by (−1/9)G I=1 (t) [41].
The I = 1 finite-size corrections are given in Table VII for every ensemble. They are computed as in [20], assuming a GS parametrization of the pion form factor. However, in contrast to [20], the parameters of the GS parametrization are obtained either by fitting the tail of the TMR correlator using the relations between the (E n , Z n ) and the pion form factor [18,19], or by using the results for m ρ and g ρππ from a dedicated pion form factor calculation, when available. This concerns ensembles C101, N401, N203, N200, D200 and J303.
We have neglected finite-size effects for the connected strange contribution, except for the SU(3) symmetric ensembles, where finite-size effects are the same as for the light-connected contribution 7 . Similarly, no finite-volume correction is applied to the charm-quark contribution.
We have performed a direct lattice calculation of the FSE on two ensembles, N101 and H105, with different volumes, L = 2.8 fm and 4.1 fm, at a common pion mass of 280 MeV. Figure 5 shows that a finite-size effect is clearly visible and statistically significant. After the finite-size correction obtained via the GS model for the pion form factor, the two correlators are in excellent agreement. This test gives us confidence that the finite-size correction we apply is reliable at our level of statistical precision.

D. The quark-disconnected contribution
We have computed the quark-disconnected contribution on a number of lattice ensembles, namely H105, N401, N203, N200, D200, N302. A typical integrand is shown in Fig. 6. The signal for the quark disconnected contribution is lost around t = 1.5 fm. Given that the absolute error of the integrand for a hvp µ grows asymptotically, it is clear that additional information constraining the tail of the disconnected TMR correlator is mandatory.
We have therefore adopted the following strategy. In our N f = 2 + 1 simulations, the isoscalar correlator G I=0,c / (t) of the (u, d, s) quarks 8 admits a positive spectral representation analogous to Eq. (13), with positive prefactors multiplying the exponentials. We expect that on the ensembles on which we have computed the disconnected diagrams, the dominant exponential in a large window of Euclidean times corresponds to the ω meson mass. As we did not perform a dedicated calculation of the ω mass, we use our determination of the ρ resonance mass. Since the latter is slightly lower than the ω mass, this is a conservative choice. We can therefore apply the bounding method in the following form, In order to quote a value a hvp,disc µ for the quark-disconnected contribution to a hvp µ , we subtract the connected light and strange contributions from the isoscalar contribution a Our results for a hvp,disc µ are listed in Table V in Appendix A.

IV. RESULTS AT THE PHYSICAL POINT
Having determined the various contributions to a hvp µ on a number of gauge ensembles, we proceed to extrapolate these results to the continuum and to the physical pion mass, m π = 134.97 MeV. We use as chiral expansion variable the dimensionless ratio where m π and f π have been determined on each ensemble.

A. The connected strange and charm contributions
For the strange-quark contribution, the statistical error (excluding the lattice spacing uncertainty) is below 1% for all the ensembles, and in many cases below 0.5%, typically for those ensembles with close-to-physical quark masses. See Table IV. The error is therefore dominated by the scale-setting uncertainty, which enters through the combination tm µ in the integrand (1). We extrapolate the results of the individual ensembles to the physical point using the fit ansatz a hvp,s µ (a, y, d) = a hvp,s µ (0, y exp ) + δ d a 2 + γ 1 ( y − y exp ) + γ 2 ( y log y − y exp log y exp ).
The index d labels the discretization, local-local or local-conserved. We observe a rather mild continuum extrapolation and both discretizations are in very good agreement. The fit goes perfectly through our physical mass ensemble and our final result for the connected strange-quark contribution is a hvp,s µ = (54.5 ± 2.4 ± 0.6) × 10 −10 , where the first error is statistical and the second is the systematic error from the chiral extrapolation. The latter is estimated from the difference between the results obtained if one includes or excludes ensembles with m π > 300 MeV. The chiral and continuum extrapolation is illustrated in the left panel of Fig. 7.
For the charm-quark contribution, the statistical error is below 0.3% for all the ensembles, and the error on the tuning of the charm hopping parameter is of similar magnitude. The error is again dominated by the scale-setting uncertainty. As can be seen on the right panel of Fig. 7, the lattice discretization of the correlator using two local vector currents leads to large cut-off effects: we observe a discretization effect of almost 70% at our coarsest lattice spacing. By contrast, for the local-conserved discretization the discretization effect is only 8%. Thus we prefer not to use the local-local discretization in our continuum extrapolation of the connected charm contribution. Furthermore, the data also suggest a very flat chiral behaviour, and we therefore use the fit ansatz a hvp,c µ (a, y) = a hvp,c µ (0, y exp ) + δ a 2 + γ 1 ( y − y exp ).
At the physical point, we obtain where the first error is statistical and the second is the systematic error induced by the chiral extrapolation. The chiral and continuum extrapolation is illustrated in Fig. 7 (right panel). A comparison of the strange and charm contributions to a hvp µ with recent publications is shown in Fig. 10.

B. The connected light-quark contribution
We have achieved a statistical error of just over two percent on a hvp,l µ on the physical-mass ensemble E250, and of 1.0 − 1.2% on all other ensembles. An important role of the other ensembles is to constrain the continuum limit, which would be very costly to achieve directly at the physical pion mass. Our lattice data points are displayed as a function of y in Fig. 8, with and without the rescaling of the muon mass with f π . We observe that the rescaled data on the right panel has a reduced dependence on y, as well as on the lattice spacing. We therefore decide to use the rescaled data for our primary analysis, but also perform the analysis of the unrescaled data in parallel for comparison.
The expected chiral behaviour of the light connected contribution is reviewed in section II E. Taking into account these considerations, we have used the following ansätze to simultaneously extrapolate our results to the continuum and to physical quark masses: where d is a label for the local-local or local-conserved correlator. All ansätze contain four parameters to be fitted, including an O(a 2 ) term to account for discretization errors. Ansatz (b) assumes a purely polynomial behaviour in the variable y, while fit (d) allows for a nonanalytic y log y term. The latter ansatz was used in our previous N f = 2 calculation [20]. Ansätze (a) and (c) are directly motivated by the discussion in section II E, (a) containing the logarithmic singularity that appears in the limit m π → 0 at fixed muon mass, while (c) contains the 1/m 2 π term relevant in the regime m µ m π m ρ . We give the results we obtain from these four ansätze, with and without rescaling m µ , in Table III. We have performed these fits either including all ensembles, or imposing cuts on y, corresponding to pion masses below 360 MeV or, alternatively, below 300 MeV. Focusing first on the rescaled data, we note that fits (a), (c) and (d) yield χ 2 /d.o.f. ≈ 1.0 while fit (b) produces higher values of around 1.6. With the pion-mass cut at 360 MeV, one sees that results (a) and (c) show good consistency and yield somewhat larger values of a hvp,l µ than fits (b) and (d). Given the more singular chiral behaviour of ansätze (a) and (c), this outcome is not unexpected. Looking at the stability of the final value for a hvp,l µ as a function of the pion-mass cut, we observe excellent stability in the case of fits (a) and (c), while the results of fits (b) and (d) systematically drift upward as a stronger pion-mass cut is imposed. With the strongest cut, m π < 300 MeV, all four ansätze yield the same result within half a standard deviation. In view of the greater stability of fits (a) and (c) against pion-mass cuts, and the stronger theoretical motivation underlying them, we choose to average the results of fit (a) and (c) with the cut m π < 360 MeV for our final central value. As a systematic error, we take the full difference between the results of these fits, and thus our final result for the connected light-quark contribution is a hvp,l µ = (674 ± 12 ± 5) × 10 −10 .
A few further remarks are in order. It is important to note that the results of fits (a) and (c) are in very good agreement with the values of a hvp,l µ directly obtained on ensemble E250 with the rescaled muon mass; see Table VI. We also remark that the statistical uncertainty on the final result Eq. (29) is only 20% lower than the statistical uncertainties on E250; we conclude that the chiral extrapolation of our results obtained at heavier pion masses, which tend to be more precise, does not lead to an artificially small final uncertainty. A comparison with the extrapolated results obtained from the standard kernel, shown in the left part of Table III, shows that the latter lie systematically higher than the rescaled ones.

C. The quark-disconnected contribution
The quark-disconnected contributions have been computed on a subset of the gauge ensembles, as described in Section II B. Three ensembles at the same lattice spacing -N203, N200 and D200 -allow us to study the chiral behaviour. Two other ensembles, N401 and N302, enable us to constrain the discretization effects.
The quark-disconnected contribution vanishes exactly for the ensembles generated at the SU(3) symmetric point. In fact, it is a double zero in the SU(3) breaking combination (m s − m l ). Since our ensembles follow a chiral trajectory at fixed bare average quark mass (2m q,l + m q,s ), we can consider the values of a disc µ as being to a good approximation 9 a function of the single variable m 2 K − m 2 π . The results of all five ensembles are thus displayed in Fig. 9 as a function of ∆ 2 2 , where ∆ 2 ≡ m 2 K −m 2 π , since close to the SU(3) symmetric point, the dependence of a disc µ on ∆ 2 2 is linear. We observe that discretization effects are negligible at the current level of precision. The result of an extrapolation to the physical point ∆ 2 = 0.227 GeV 2 assuming a linear proportionality to ∆ 2 2 is a hvp, disc µ = −18.6(2.2) × 10 −10 . As discussed below Eq. (18), the disconnected contribution has a singular behaviour in the limit m π → 0, closely related to the corresponding behaviour of the connected light contribution. Therefore, we consider the possibility that the disconnected contribution contains a term with precisely the dependence given in Eq. (18). In order to make this term consistent with the double zero of the disconnected correlator at ∆ 2 = 0, we fixM 2 ≡ 1 2 m 2 π + m 2 K to its physical value, express m 2 π through the variable ∆ 2 and use the ansatz Fitting the single free parameter γ 8 , we obtain a hvp, disc µ = −27.7(2.2) × 10 −10 . From Fig. 18, it is clear that both the linear fit in ∆ 2 2 and the one based on ansatz (30) are consistent with the lattice data. While a singular chiral behaviour must be present in a hvp, disc µ , the ansatz (30) may lead to an overestimate of this effect. Therefore, we quote as our final result the average of the linear and the chirally singular fit, where the first error is statistical and the second is a systematic error associated with the extrapolation to the physical point, taken to be the half-distance between the two extrapolated values.
where the first error is statistical and the second is the systematic error. The latter is dominated by the chiral extrapolation of the light-connected and the disconnected contributions. The result Eq. (32) does not contain any correction for QED or strong isospin-breaking effects. For now, we do not attempt to include such a correction, but rather add (in quadrature) a systematic uncertainty of 7.2 × 10 −10 corresponding to a recent lattice calculation of these effects [54]. This then leads to our final result given in Eq. (33) below. 9 A residual dependence on the independent combination ( 1 2 m 2 π + m 2 K ) persists at higher orders in the chiral expansion and via O(a) discretization effects. In the rightmost panel, the full results, including (where available) the contributions from quark-disconnected diagrams and corrections due to isospin-breaking, are compared to the phenomenological determination of Ref. [55], represented by the red vertical band. Our result is compared to the calculations labelled FNAL-HPQCD-MILC 19 [56][57][58], PACS 19 [59], ETMC 19 [54,60,61], RBC/UKQCD 18 [39], BMW 17 [38], as well as our previous calculation in two-flavour QCD [20] (Mainz/CLS 17).

V. DISCUSSION AND COMPARISON
In this paper we have presented a calculation of the hadronic vacuum polarization contribution to a µ based on gauge ensembles with N f = 2 + 1 flavours of O(a) improved Wilson quarks. Our final result is a hvp µ = (720.0 ± 12.4 stat ± 9.9 syst ) · 10 −10 , where the first error is statistical, and the second is an estimate of the total systematic uncertainty, which also accounts for the fact that the corrections due to isospin breaking have not been included. We thus find that the overall error of our determination is 2.2%. In Fig. 10 we compare our results to those of several other recent lattice calculations [20,38,39,54,58,59]. While our estimate is at the higher end of lattice results, we note that the direct difference with the result based on dispersion theory of Ref. [55] is 26.6 ± 16.0, which amounts to ∼ 1.7 standard deviations and may signal a slight tension. There are several ways in which our result can be improved without relying on the obvious strategy of adding more ensembles and increasing the overall statistics. First, we have seen in section III B that the use of detailed spectroscopy information in the isovector channel is a huge advantage, as it nearly halves the statistical uncertainty in the estimate for a hvp,l µ on ensemble D200. This is the result of either constructing the vector correlator from the energies and overlaps determined via the GEVP or of using this information in the improved bounding method. Extending these calculations to more ensembles -in particular those with physical and near-physical pion masses -will boost the statistical accuracy and reliability significantly.
Second, we have pointed out that it is advantageous to split the correlator into isovector and isoscalar components according to Eq. (4) rather than focussing on separating the contributions from individual quark flavours. One reason is that the singular chiral behaviour expected from Eq. (17) is shared between the light quark connected and the disconnected contributions. This will help to better constrain the pion mass dependence of the quarkdisconnected contribution, which is often still obtained from an extrapolation to the physical point from a set of results at heavier pion masses. The decomposition according to isospin also gives a better handle on finite-volume effects, which are partly compensated between the light connected and disconnected contributions. This is of particular importance, since finite-volume corrections for the disconnected part of the vector correlator could be sizeable but, to our knowledge, have not been estimated so far.
The third refinement concerns the determination of isospin-breaking corrections. We stress again that our final estimate in Eq. (33) is valid at a well-defined reference point of the isospin-symmetric theory, given by the mass of the neutral pion in the continuum limit. The determination of the corrections due to isospin breaking relies on the definition of an alternative reference point that is consistent with the effects induced by a non-vanishing mass splitting among the up and down quarks and the coupling between quarks and photons. This requires an adjustment of bare parameters and the re-evaluation of a number of observables that enter the calculation of a hvp µ . An account of the status of our activities in this direction is given in Refs. [14,15]. In the absence of a complete evaluation, we have refrained from simply adding results for the isospin-breaking correction from the literature. Instead, we have opted for an additional systematic error which is as large as the correction determined in [54].
As the community awaits the first results from the E989 experiment at Fermilab, it is remarkable that several collaborations using different setups and discretizations of the QCD action obtain largely consistent estimates for a hvp µ with overall errors at the level of 2%. However, the collection of available results does not allow for a firm conclusion as to whether the phenomenological estimate or the so-called "No New Physics" scenario is confirmed.
ensembles.  Values of the connected light contribution to a hvp µ , in units of 10 −10 , for the local-local (LL) and for the local-conserved (CL) discretizations of the correlation function, as described in the main text. The symmetric lattice derivative of the improvement term is used. In addition, the pion mass and the pion decay constant are given, some of them taken from Ref. [64]. No FSE correction has been applied on any data in the table. The treatment of the long-distance part of the correlator is described in section III A. In the "No rescaling" columns, the scale-setting error has not been included. In the "With rescaling" columns, the statistical fluctuations of the f π determination are taken into account.  ≡ a hvp µ (L = ∞) − a hvp µ (L) = ∆a < µ (t i ) + ∆a > µ (t i ) on the isovector contribution a hvp,I=1 µ in the TMR in units of 10 −10 . We used the value t i = (m π L/4) 2 /m π , as in [20], to which we refer for unexplained notation. The parameters of the GS model are obtained by fitting the tail of the correlation function using the Lüscher formalism. When there are two lines, the second line corresponds to the GS parameters extracted from a direct lattice calculation of the timelike pion form factor [29]. The FSE effects are given both with and without rescaling of the muon mass via f π .

GS parameters
No rescaling With rescaling id We list in Table IX the values of the cs pseudoscalar meson masses determined for different values of the charm hopping parameter. The 'physical' value of the charm-quark hopping parameter is determined by the condition that the cs meson mass match the physical value of the D s meson mass, m Ds = 1972 MeV.  (19) and (20). For ensemble J303, the data from [29] has been rebinned and analyzed using the jackknife method.