Hadronic vacuum polarization contribution to the anomalous magnetic moments of leptons from first principles

We compute the leading, strong-interaction contribution to the anomalous magnetic moment of the electron, muon and tau using lattice quantum chromodynamics (QCD) simulations. Calculations include the effects of $u$, $d$, $s$ and $c$ quarks and are performed directly at the physical values of the quark masses and in volumes of linear extent larger than $6\,\mathrm{fm}$. All connected and disconnected Wick contractions are calculated. Continuum limits are carried out using six lattice spacings. We obtain $a_e^\mathrm{LO-HVP}=189.3(2.6)(5.6)\times 10^{-14}$, $a_\mu^\mathrm{LO-HVP}=711.1(7.5)(17.4)\times 10^{-10}$ and $a_\tau^\mathrm{LO-HVP}=341.0(0.8)(3.2)\times 10^{-8}$, where the first error is statistical and the second is systematic.

We compute the leading, strong-interaction contribution to the anomalous magnetic moment of the electron, muon and tau using lattice quantum chromodynamics (QCD) simulations. Calculations include the effects of u, d, s and c quarks and are performed directly at the physical values of the quark masses and in volumes of linear extent larger than 6 fm. All connected and disconnected Wick contractions are calculated. Continuum limits are carried out using six lattice spacings. We obtain a LO-HVP Introduction.-Ever since the discovery of the electron's spin [1,2], the magnetic moments of leptons have accompanied the development of quantum mechanics and quantum field theory. This is particularly true of the small, "anomalous," quantum corrections to these moments, a , where denotes either the electron (e), the muon (µ) or the tau (τ ) (see e.g. [3] for an introduction). Today, a e is one of the most precisely measured [4] and computed [5,6] quantities in nature, with a total uncertainty below 1 ppb. Theory and experiment agree and the measurement can be used to make the most precise determination of the fine-structure constant α [6].
In the case of the muon, the precision of the measurement [7] and of the standard model (SM) prediction (e.g. [8]) are roughly matched at around 0.5 ppm. However, theory and experiment disagree by more than 3 standard deviations. This is particularly enticing, because it could be a sign of new, fundamental physics. The anomalous magnetic moment of the muon is generically much more sensitive to new, massive degrees of freedom than that of the electron. This is because, in many extensions of the SM, the contributions of new particles are proportional to the lepton mass squared, which is roughly 4 × 10 4 times larger for the muon. Moreover, a new experiment is beginning to take data at Fermilab [9], with the goal of reducing errors by a factor of four, and another one is planned at J-PARC, with similar objectives [10].
The same argument should make a τ even more interesting for new physics searches: the τ mass is close to 17 times that of the muon. However, its very short lifetime, of order 10 −13 s, has meant that no direct measurement of a τ has yet been made, though a concrete proposal for doing so [11] is being implemented [12].
Theoretically, the leading source of uncertainty in the SM prediction of a µ is the leading order (LO) hadronic vacuum polarization (HVP) contribution, a LO-HVP µ , which is responsible for over 79% of the total error [8].
Here we present lattice QCD calculations of the LO-HVP contribution to the anomalous magnetic moments of all three leptons. The calculations include all contributions from u, d, s and c quarks, directly at the physical values of their masses, in their quark-connected and quark-disconnected configurations. Contributions from third generation quarks can easily be estimated and are found to be much smaller than our statistical errors, even for the τ that is most sensitive to them (see e.g. [40] for a calculation of the b contribution to a LO-HVP µ ). Some previous lattice calculations of a LO-HVP , =e, µ [25,31], and of a LO-HVP , =e, µ, τ [41], included all of these contributions, but involved difficult extrapolations to the physical value of the average u-d quark mass and only estimated the disconnected parts. In the present Letter, we work directly at the physical point and compute disconnected contributions directly. Moreover, we implement a description of the lattice results [42,43] that solves the small virtuality issue [44] with finite-volume (FV) artifacts that are exponentially suppressed in lattice size.
A unified treatment of the HVP contribution to the three lepton anomalous magnetic moments provides important cross-checks that validate the methods used. As the typical virtualities probed by these moments are around m 2 /4, the vast difference in the mass of the leptons means that a large range of relevant scales are checked. In particular, agreement between our results and phenomenology in the electron case validates our understanding of small virtualities, and of larger virtualities in the τ case. In addition, the inclusion of all flavors up to the charm allows a controlled matching onto perturbation theory. Thus, all energy scales from zero to infinity are controlled in our calculation.
Methodology.-We consider the zero three-momentum, two-point function of the quark electromagnetic current in Euclidean time t: with e the positron charge, x=(t, x) and j µ /e= 2 3ū γ µ u − 1 3d γ µ d − 1 3s γ µ s + 2 3c γ µ c. We work in the isospin limit, m u =m d . Because C µν 's flavor components are calculated separately and have different statistical and systematic uncertainties, it is useful to treat them separately. Physically, an isospin separation is useful. Thus, where in the top equality the first three terms correspond to the quark-connected contractions of the light (u and d combined), strange and charm quarks, and the fourth to the quark-disconnected contractions of all four flavors.
In the second equality, the separation is made between isospin I=1 and I=0 contributions, given by C I=1 µν = 9 10 C ud µν and C I=0 µν = 1 10 C ud µν + C s µν + C c µν + C disc µν . It is straightforward to obtain the corresponding LO-HVP contributions to the anomalous magnetic moment of lepton from these correlation functions [21,45,46]: with ω(r) = π 2 r + 2 − r(r + 4) 2 / r(r + 4), α = e 2 /(4π) and where the scalar polarization function renormalized in the Thomson limit is given by (see also [42]) In Eqs. (3) and (4), the superscript f can stand for ud, s, c, disc, I=1, I=0 and where the " " indicates that this equation also applies to the full LO-HVP contribution. Eq. (4) implicitly includes the subtraction of the polarization tensor Π µν (Q=0), which was shown in [28] to be critical for reducing FV effects and, through the factor t 2 , the subtraction of the polarization scalar Π f (0). On a T × L 3 lattice with spacing a, the integral over t in Eq. (4) is replaced by a sum, in increments of a, that runs up to T /2, once the correlator C f ii (t) has been averaged with C f ii (T − t). Moreover, the integral over Q in Eq. (3) should, in principle, be replaced by a sum from 0 to π/a in steps of 2π/T . Here we keep the integral, but cut it off at a value Q=Q max , chosen much smaller than π/a, so as to keep discretization errors under control, but above which perturbation theory can be applied. Then we decompose the anomalous magnetic moments of the leptons into three terms: where the low momentum contribution, a LO-HVP ,f (Q≤Q max ), is obtained from the lattice as described above, and where the last term is the highmomentum, contribution renormalized at Q max and computed in perturbation theory [47]. The second term in Eq. (5) is required to shift the renormalization point from Q max to Q=0. It is obtained with lattice results for C f ii (t), through Eq. (4) with Q=Q max . γ (Q max ) is a known kinematical factor [47]. In obtaining Eq. (5), it is assumed that ∆ pert a LO-HVP ,f (Q>Q max ) is equal to the value that it would have nonperturbatively. We check this by studying the dependence of our results on the choice of Q max .
The replacement of the FV sum over Q by the corresponding integral is our choice of interpolation for the HVP functionΠ(Q 2 ). It constitutes an alternative to e.g. the Padé approximation proposed in [44], the Marichev interpolation advocated in [48] or the finite-energy sum rule approach of [49]. The integrand in Eq. (3) has no singularities in the region of integration. Therefore, an application of Poisson's summation theorem guarantees that the corrections entailed in replacing the sum by an integral over Q are exponentially small in T . These corrections will be accounted for in our estimate of FV uncertainties.
Lattice details.-We employ a tree-level improved Symanzik gauge action [50] and a fermion action for N f =2 + 1 + 1 flavors of stout-smeared [51], rooted, staggered quarks. We have generated 15 ensembles at six values of the bare coupling, β, corresponding to lattice spacings ranging from 0.064 to 0.134 fm. The average up and down quark mass and the strange quark mass are tuned to around the physical mass point defined using the Goldstone pion and kaon masses. The charm quark mass is fixed in units of the strange mass to m c /m s =11.85 [52]. The spatial dimensions of our lattices are in the range 6.1-6.6 fm and the temporal ones in the interval 8.6-11.3 fm. The lattice spacing is fixed with the pion leptonic decay constant, f π . At each value of β, between 450 and 3500 configurations, separated by 10 unit length rational hybrid Monte Carlo (RHMC) [53] trajectories, are used. Details are given in [47] and more information about the simulations can be found in [54].
For the electromagnetic current correlator, we use the conserved lattice current at the source and sink so that no renormalization is necessary. We calculate the connected contributions to the correlators using point sources. We use the all-mode-averaging (AMA) technique of [55] and 768 random source positions on each configuration for the light quarks, 64 sources for the strange and 4 for the charm. To compute the quark-disconnected contributions, we apply AMA again, and exploit the approximate SU(3) flavor symmetry on around 6000 stochastic sources [30,56]. These are random, four-volume sources with which we compute the zero-momentum, time propagators, correcting for bias. For the disconnected contribution of the charm we apply a hopping parameter expansion.
Analysis.-Even with our high statistics, the signal deteriorates quickly with increasing distance in our light and disconnected correlators. Thus, in implementing Eqs. (3) and (4), we introduce cuts, t c , in time beyond which we replace the correlator by the average of an upper and a lower bound [35,47]. t c is chosen such that the upper and lower bounds agree well within statistical errors and where these errors are not too large. The upshot is that our result for the light contribution to a LO-HVP is obtained by summing the integrand in Eq. (4) with C ud ii (t) given by our lattice data up to t c and performing the rest of the sum from t>t c to T /2 with C ud ii (t) replaced by the bound average. The results of this procedure for t c in the range of (3.000±0.134) fm are averaged to account for possible statistical fluctuations in the correlator at a given t c . The disconnected contribution to a LO-HVP is obtained in an identical fashion, but with t c in the range (2.600 ± 0.134) fm.
We limit the integral over Q in Eq. (3) to Q max , and use perturbation theory to obtain the complement. We consider Q 2 max =1, 2, · · · , 5 GeV 2 . In what follows, quantities with the subscript "lat" correspond to lattice results obtained in a given simulation. Their dependence on lattice spacing and quark masses will be left implicit. To extrapolate our results a LO-HVP ,f,lat (Q≤Q max ) to the continuum limit and to interpolate them to the physical mass point, we fit them to a function which depends on the Goldstone pion and kaon masses squared, on the η c mass and on the lattice spacing squared [47]. Since the simulations are performed close to the physical mass point, a constant or linear dependence in the mass parameters is always sufficient. Moreover, for all flavor contributions, good fit qualities can be achieved with a linear a 2 dependence for all three leptons and all values of Q max considered here. Because taste violations play an important role in the continuum extrapolation of a LO-HVP ,ud , we have also tried correcting for these effects using one-loop staggered chiral perturbation theory before performing a continuum extrapolation [34]. While the continuum extrapolation is significantly milder, the continuum limit results obtained are consistent with the ones presented here. Our continuum extrapolations are discussed in detail and examples are shown in [47]. Here we emphasize that with simulations at six lattice spacings down to 0.064 fm, we have full control over the continuum extrapolations.
This analysis yields the continuum extrapolated flavor quantities, a LO-HVP ,f (Q≤Q max ), for the five values of Q max considered. For each value of Q max , we sum the appropriate flavor quantities, to get the corresponding I=1, I=0 and total, low-Q contributions to the lepton anomalous magnetic moments. The results for these and the individual flavor contributions are given in [47], with statistical and systematic errors obtained as described below. To these contributions we add the corresponding complements given in Eq. (5). These complements require the computation ofΠ f (Q 2 max ). This is done using Eq. (4) and requires a continuum limit and physical mass point interpolation very similar to that performed for a LO-HVP ,f,lat (Q≤Q max ) [47]. The perturbative contributions, ∆ pert a LO-HVP ,f (Q>Q max ), are computed from results for Π f (Q 2 ), with terms up to O(α 4 s ), obtained using the code rhad [57], as explained in [47]. These corrections are below our statistical errors for the e and µ, which have very little sensitivity to large Q, but are significant for the τ . In [47] we study the Q max dependence of our results for a LO-HVP ,f . The fact that they are independent of Q max ≥ √ 2 GeV within errors, in particular for =τ , indicates that our continuum-limit, lattice results are consistent with five-loop perturbation theory for momenta, Q, above that value.
Systematic errors and results.-The procedure described above yields a LO-HVP ,f for all f and for all three leptons. In our physical fits, the errors associated with the small interpolations in mass are negligible. Those associated with the continuum extrapolations are not. To estimate them, we impose four cuts on the lattice spacing: no cut, and a≤0.118, 0.111, 0.095 fm. This number is reduced to three in the disconnected case for which we have no results at a = 0.064 fm. The systematic error associated with the matching to perturbation theory is determined from the results with Q 2 max =2, · · · , 5 GeV 2 , covering the range safely accessible to our lattice calculations and to perturbation theory. The one associated with the time cut is determined by considering t c ranges shifted by −0.134 fm compared to those given above. Over the total range of t c considered, the two-pion bounds change by a factor close to 3. The final central value is the unweighted average of all results. Each systematic error component is chosen to cover all of the central values resulting from the variation, over the ranges described above, of the variable associated with that component. Furthermore, we add a 0.8% systematic error to our results for a LO-HVP ,f (Q≤Q max ) due to the uncertainty in our determination of the lattice spacing [47]. The statistical error is the jackknife error of the central value over jackknife samples with bins of length 10 configurations. The results for the individual flavor contributions to the magnetic moments of all three leptons are given in [47].
In the absence of a systematic study with simulations in a variety of volumes, only model estimates of FV effects can be made. As argued in [58,59], for large volumes those effects will be governed by pion contributions that can be computed in chiral perturbation theory (χPT) [58]. Since the I=0 channel is dominated by = e (units of 10 −14 ) = µ (units of 10 −10 ) = τ (units of 10  LO-HVP contribution to the anomalous magnetic moments of the e, µ and τ leptons. The first two lines give our results for the I=1 and I=0 contributions. The I=1 results include the FV corrections, which are negligible in the I=0 case. The last line displays our results for the total LO-HVP contribution. In addition to the terms included in the I=1 and I=0 components, this total also accounts for QED and m d =mu corrections. The first error on all results is statistical, the second is associated with the continuum extrapolation, the third with our bounding procedure, the fourth with the matching to perturbation theory, the fifth with the the lattice spacing uncertainty and, where applicable, the sixth with the FV correction and the seventh with the IB correction. three-pion exchange, the FV effects are expected to be smaller than those of the I=1 contribution, which are already small. Thus we consider only the latter. Our computation of these effects is summarized in [47] and the appropriate corrections are added to our I=1 and total results. They are 4.6(4.6) × 10 −14 for the e and 13.5(13.5) × 10 −10 for the µ with negligible Q max dependence in the range of interest. For the τ they range from 9.4(9.4)×10 −9 to 1.6(1.6)×10 −8 for Q max = 1÷ √ 5 GeV. We associate with these corrections a 100% uncertainty included in our error budget.
We quote our final results for a LO-HVP for all three leptons in Table I. Combining all errors in quadrature, we obtain a LO-HVP e with an uncertainty of 3.3%, a LO-HVP µ of 2.7% and a LO-HVP τ of 1.0%. Not surprisingly, the relative error increases with the sensitivity of the anomalous moment to long-distance physics.
Discussion.-It is interesting to compare these results with those in the literature. There are only two lattice QCD determinations of the LO-HVP contribution to the muon anomalous moment which include the contributions of quarks up to the charm [25,34]. Compared to those, our calculation is the only one in which the continuum extrapolation is performed directly at the physical mass point and which includes a reliable determination of the quark-disconnected contribution. There exist also many precise phenomenological determinations of a LO-HVP µ , as discussed in the Introduction. Here we consider three recent ones [8,19,20].
We plot all of these results in Fig. 1, together with ours. Also shown on this plot is the value that a LO-HVP µ would have to have to explain the experimental measurement of a µ [7], assuming that all other SM contributions are unchanged, i.e. assuming no new physics (NP). Using the SM contributions summarized in [8], we find a LO-HVP µ,noNP = (720.0 ± 6.8) × 10 −10 . The errors on the lattice results, which are in the range of 2.0% to 4.1% are substantially larger than those of the phenomeno- with the only other two N f = 2+1+1 lattice QCD calculations [25,34] and with recent ones obtained from phenomenology [8,19,20]. For the lattice results, the first error is statistical and the second is the total error, including systematics. The shaded region is the value that a LO-HVP µ would have to have to explain the experimental measurement of aµ, assuming no new physics. logical approach. Our result for a LO-HVP µ is larger than those of the other lattice calculations and in slight tension with the one from HPQCD [34] which is 1.9 σ away. A more detailed flavor-by-flavor comparison is given in [47]. However, our result is consistent with those from phenomenology within about 1 standard deviation, as well as with a LO-HVP µ,noNP . Thus, one will have to wait for the next generation of lattice QCD calculations to confirm or infirm the larger than 3 σ deviation between the measurement of a µ and the prediction of the SM based on phenomenology.
Regarding a LO-HVP e , there are two other lattice calculations [34,41] and only one concerning a LO-HVP  [13]. Comparing these results to ours in Table I, we find the following. The result of [34] for a LO-HVP e displays a tension with ours which is slightly smaller than the one for a LO-HVP µ . On the other hand, our results are fully compatible with the phenomenological ones, indicating that we control the physics of the HVP over full range of Q 2 . In addition, our result for a LO-HVP e has an error which is about half that of the lattice result of [41] and for a LO-HVP τ , it is approximately 3 times more precise. In fact, our result for the latter is more precise than the phenomenological one.
Acknowledgments Note added.-After submission of this Letter for publication, a preprint reporting on an N f = 2+1 lattice QCD calculation of a LO-HVP µ appeared [60]. That calculation includes a lattice computation of many isospin breaking effects. Its result for a LO-HVP µ is in excellent agreement with ours.   (1) FIG. S1. Positions of our 15 simulations in the m ud -ms plane, as represented by (Mπ/fπ) 2 and (M 2 K − M 2 π /2)/f 2 π , respectively, normalized by their physical values. Here Mπ and MK are the masses of the Goldstone taste partners. At a given value of β, the polygons delimited by the points (M 2 π , M 2 K − M 2 π /2)/f 2 π indicate whether these points bracket the physical point. The numbers in parentheses, next to the values of β, indicate the number of independent simulations at that β.

Lattice parameters
The work presented here is based on 15, N f = 2 + 1 + 1 simulations performed around the physical mass point and at six values of the bare coupling, on lattices of spatial extents larger than 6 fm. The parameters of these simulations are summarized in Table S1, for each of the six bare coupling values. Because we are using staggered fermions, our results suffer from taste artefacts. The corresponding splittings in the pion spectrum are also given in Table S1.
All of our simulations are performed around the physical mass point in the isospin limit, with m d = m u and α = 0. In Fig. S1 we show the positions of our 15 simulations in the m ud -m s plane, as represented by (M π /f π ) 2 and (M 2 K − M 2 π /2)/f 2 π , respectively. Here M π and M K are the masses of the Goldstone taste partners. For four of our six lattice spacings, we have three or more different values of m ud and m s . As the figure shows, in most cases these simulations bracket the physical point and the values of (M π /f π ) 2 and (M 2 K − M 2 π /2)/f 2 π are never more than 3, respectively 6, percent away from their physical values.
To obtain a precise signal for the connected-light and quark-disconnected contributions to the current-current correlator at the large distances required for the com-putation of the light lepton magnetic moments, a significant statistics is needed. Thus, we have generated a large number of configurations, separated by 10 unit length RHMC trajectories, and have performed the computations using a large number of sources on each configuration. The number of sources required for the strange contribution is significantly smaller, because of the better quality of its signal. All of these numbers are also specified in Table S1. Moreover, as mentioned in the main text, the precision of the connected-charm contribution is so high, that we have chosen to evaluate it on a reduced set of configurations, to bring its statistical error more in line with those of the other contributions and to allow for fits with good fit qualities. We have found that 40 configurations, maximally separated in molecular dynamics time, and 4 sources are sufficient.

Time cuts and bounds for the light and disconnected contributions
In the case of the light and disconnected correlators, C ud ii (t) and C disc ii (t), i = 1, 2, 3, the signal deteriorates quickly with increasing time. Thus, we introduce a cut t c in time, beyond which we replace the correlator by an upper and a lower bound [35]. In the continuum limit, C ud ii (t) is a sum of falling exponentials with positive coefficients. Therefore, it is bounded below by zero. Moreover, because it is proportional to the isospin triplet correlator, whose lowest-energy contribution comes from a two-pion state, it is bounded above by that contribution. Therefore, in the continuum limit and the limit of infinite statistics, the correlator satisfies and E 2π is the energy of two pions, each with the smallest nonvanishing lattice momentum, for which we use 2π/L. Eq. (S1) holds up to higher-order, wrap-around contributions which are suppressed exponentially in nonvanishing multiples of T E 2π . In Fig. S2, for an ensemble at β = 3.9200 and for all three leptons, we plot the resulting upper and lower bounds on a LO-HVP ,ud (Q≤2 GeV) as a function of t c . These are obtained by inserting the bounds of Eq. (S1) into Eqs. (3) and (4) from the main text, for t>t c . For t≤t c , it is the correlator C ud ii (t) obtained directly in the simulation which is used. As the figure shows, the upper and lower bounds typically agree for t c β a [fm] (T ×L/a 2 ) a 2 ∆taste #conf-ud/#conf-s/#conf-c/#conf-disc #src-ud/#src-s/ TABLE S1. Parameters of the 15 simulations performed. For a given bare coupling, β, the simulations differ only by a slightly distinct choice of quark masses around their physical value. Here we list the β, the lattice spacings, a, the sizes T × L, the taste splittings, a 2 ∆taste = (aMi5) 2 − (aMπ) 2 (with M5i the mass of the ξ5i taste partner and Mπ, the mass of the Golstone), the number of configurations used and the number of sources per configuration used for each flavor contribution. Note that for the ud and s contributions, all 15 simulations are used while for the charmed one, 2 simulations for β = 3.9200 are omitted and, in addition, the one at β = 4.0126 for the quark-disconnected one.
The reason for averaging over t c is to dampen the effect of possible statistical fluctuations in the value of the correlator from from one t c to the next. This is particularly useful with staggered fermions, for which even and odd times are mostly uncorrelated. Thus, for each lattice spacing, we round the number of time slices in these t c intervals to the nearest even integer. For β ≥ 3.9200, this leads to four time slices in each t c interval and to a shift of two time slices between intervals. For β < 3.9200, these figures are two and one. While this average is not strictly necessary, because the bound means between two neighboring t c 's in the relevant interval agree statistically, it does reduce overall the fluctuations in central values, leading to better χ 2 /dof, e.g. in the continuum extrapolations.
The disconnected contribution alone can be constrained for large enough time separations, where the isospin singlet channel, dominated by three-pion states, can be neglected compared to the triplet one, dominated by two-pion states. Here we have, for large enough t c , Pion-pion interactions change the smallest two-pion momentum from 2π/L in that channel. Using the model of [61] and neglecting four-pion contributions, we determine the change in the two-pion energy to be around 2%. We checked that such a reduction of the momentum changes the result on a LO-HVP   (27) GeV of [63]. To convert the lattice results into physical units, we use the pion decay constant obtained from pion leptonic decays, which is free of electromagnetic corrections and, to very good accuracy, equals the decay constant in the m d = m u limit [64]. This yields a well defined physical point in the isospin limit. In intermediate steps of the analysis, we use the Wilson-flow-based [65] w 0 -scale [66] as described below.
As noted in [42,67], and generalized here for individual flavor contributions, Thus, for very light leptons, the leading dependence of a LO-HVP ,f on lattice spacing is quadratic. In practice, we find this to be almost exactly true for the electron. For the µ, this dependence ranges from quadratic to about 20% weaker than quadratic, depending on the flavor contribution. Only for the τ can it be significantly less for some contributions, as low as a 0.8 . Because it maximizes the error bar due to the lattice spacing uncertainty, we will assume a quadratic dependence on a for all leptons and for all flavor contributions. This means that, for a given simulation, the relative error on a LO-HVP ,f,lat (Q≤Q max ), due to the scale uncertainty, will be taken to be twice that on the lattice spacing. While this error may be modified some by the combined mass interpolations and continuum extrapolations, it represents the dominant effect on our results for a LO-HVP ,f (Q≤Q max ). This strong dependence on a means that uncertainties on the lattice spacing must be estimated with care. As noted above, we fix a with f π , through a/w 0 and w 0 f π . a/w 0 is obtained to sub-permil precision for each of our simulations. On a subset of these simulations, we determine w 0 f π in the continuum limit with a total 0.4% error. Using f π = 130.50(1)(3)(13) MeV from the PDG [68,69], this means that we have a 0.4% error on our lattice spacings [70]. As discussed above, this translates into a 0.8% uncertainty on our results for a LO-HVP ,f (Q≤Q max ), =e, µ, τ , f =ud, s, c, disc, which we carry over to our final results for a LO-HVP ,f , f =(I=1), (I=0), . It is important to note that this relative error on the individual flavor contributions to a LO-HVP is 100% correlated across these contributions. Thus, it is the same for the individual contributions as it is for their sum. To obtain numbers that can be compared with experiment, we must interpolate our results to the physical mass point and extrapolate them to the continuum limit. We do so by fitting them to a function which depends on the Goldstone pion and kaon masses squared, on the η c mass and on the lattice spacing squared. Since the simulations are performed close to the physical mass point, a constant or linear dependence in the mass parameters is always sufficient. Moreover, for all flavor contributions, good fit qualities can be achieved with a linear a 2 dependence for all three leptons and all values of Q max considered here. Thus, the physical values of the contributions, a LO-HVP +γ a,f a 2 + γ π,f (M lat π /2, for = e, µ, τ and for f =ud, s, c, disc. It turns out that for the statistical accuracy reached here, the light quark mass in our simulations is sufficiently well tuned that γ π,f comes out consistent with zero for all flavors, except the charm, and omitting it increases the p-values of the fits. Thus we set it to zero for all but the charm contribution. This can be done for γ K,f in the case of the light, quark-disconnected and charm contributions. γ K,f must be kept as a free parameter for the connected strange contribution. Moreover, γ c,f is needed for the charm contribution to correct for a slight mistuning of m c . Here we focus on the results obtained with Q max =2 GeV. Projections of these fits to our simulation results for a LO-HVP all of these situations in the same way and discuss them together. For the light-quark contribution to a LO-HVP (Q≤2 GeV), the dependence on meson masses is not significant statistically and the terms associated with this dependence can be ignored. However, as can be seen in the upper panels of Fig. S4, the dependence on a 2 is strong, due to the sensitivity of this contribution to low-energy, two-pion states which, in turn, are sensitive to taste splittings. The fact that the anomalous moment of the lighter e is more sensitive to these states than that of the µ that is, in turn, more sensitive than that of the τ , explains the fact that a LO-HVP e (Q≤2 GeV) has the strongest a 2 dependence while a LO-HVP τ (Q≤2 GeV) has the weakest.
The situation is different for the strange contribution, much less affected by taste violations. As the second panels of Fig. S4 show, the continuum limits are very mild. They are much less so for the charm, as shown in the third panels, due to the large value of m c in lattice units.
Here it is the magnetic moments of the more massive leptons which are steeper, due to their sensitivity to larger Q. In addition to the dependence on a 2 , a linear dependence on M 2 K χ is needed for both contributions and one on M ηc is required to correct a slight mistuning of the charm mass in that quark's contribution.
Our results for a LO-HVP ,disc,lat (Q≤2 GeV) have large lattice artefacts, as shown in the bottom panels of Fig. S4. This is because the taste violations of the ud contribution enhance the SU(3)-flavor cancellation against the s contribution in a LO-HVP µ,disc,lat , as a 2 is increased. In these results we neglect the charm contribution which we find to be less than 1% of the total disconnected contribution on our coarsest lattice, i.e. much smaller than the disconnected, statistical error. In addition, because statistical errors are quite large, no dependence on quark mass is required to describe the lattice data.
As explained in the main text, the systematic error associated with these continuum limits and physical-point interpolations are obtained by imposing four cuts on the lattice spacing in the quark-connected case and three for the disconnected contributions. The results of these cuts are then combined, as detailed in the main text, to give a central value and statistical and systematic errors. The results for the various contributions to the magnetic moments of all three leptons with the four values of the momentum cut Q max considered are summarized in Table S2. . Results for Q 2 max = 1, · · · , 5 GeV 2 are given. In the case of the e and µ, this "low-energy" contribution is equal to the full LO-HVP contribution within errors. Only in the case of the τ is it necessary to add a "highenergy" complement at the current level of precision. In the results presented, the first error bar is statistical, the second is the systematic uncertainty associated with the continuum extrapolation and the third with the bounding procedure described in Sec. 2. The latter does not affect the s and c contributions and is denoted (−) for these contributions. Note that finite-volume corrections discussed below in Sec. 8 are not included here.

Determination ofΠ
In order to determine the various flavor f contributions to the anomalous magnetic moments of the leptons through Eq. (5), we also need to computeΠ f (Q 2 max ). Even for large Q 2 max , this is a non-perturbative quantity because it requires a subtraction at Q 2 = 0. We obtain Π f (Q 2 max ) from our lattice results for the current-current correlator C f ii (t) through Eq. (4) of the main text. We do so for each one of our simulations, for each of the flavors, f =ud, s, c, disc and for the five values of Q max considered here. The resulting valuesΠ f lat (Q 2 max ) must then be interpolated to the physical mass point and extrapolated to the continuum. To this end, we follow the exact same procedure as for a LO-HVP ,f . That is, we fit our lattice results to a functional form which parametrizes the dependence on lattice spacing and on quark masses. As for the contributions to the leptonic magnetic moments, we find that a constant or linear dependence on a 2 , M 2 π , M 2 K χ and M ηc describes the lattice results well, for all f and Q max . Again with the charm we have to reduce the statistics to get acceptable χ 2 /dof. Since the features of these fits are very similar to those for a LO-HVP ,f,lat (Q≤Q max ), we do not repeat the discussion here. Instead we plot, in Fig. S5, the results of these fits for all four flavors and for the Q 2 max =4 GeV 2 used in other plots. The results corresponding to all values of Q 2 max considered here and to all physically relevant flavor combinations are summarized in Table S3. Central values and statistical and systematic errors are obtained, as for a LO-HVP ,f (Q≤Q max ), using flat distributions, cuts in a 2 and the bounding procedure described in Sec. 2.  ), as obtained on the lattice through Eq. (4) of the main text, with a continuum extrapolation and an interpolation to the physical mass point. In our conventions,Π f (Q 2 max ) is smaller by a factor of 4π 2 than it is, for instance, in [34]. Results for different flavor combinations, f =ud, s, c, disc, I=0, , are given. Note thatΠ I=1 = 9 10Π ud . In this table, the first error is statistical, the second is the systematic uncertainty associated with the continuum extrapolation and the third with the bounding procedure described in Sec. 2. The latter does not affect the s and c contributions and is denoted (−) for those quantities. Note that finite-volume corrections discussed below in Sec. 8 are not included here.

Separation of a LO-HVP into low-and high-virtuality contributions
The computation of a LO-HVP requires a determination of the renormalized, scalar polarization function,Π(Q 2 ) for all values of the Euclidean momentum Q, from 0 to ∞. Of course at finite lattice spacing, momenta up to infinity are not available. Thus, as explained in the main text, we deal with this problem by separating the low and high virtuality contributions at a value of Q = Q max in the following way (repeating Eq. (5) of the main text): where the low momentum contribution, a LO-HVP ,f (Q≤Q max ), is obtained from the lattice as described in the main text, and where the last term is the high-momentum, perturbative contribution renormalized at Q max , i.e.
where the kinematic function, ω(r), is defined after Eq. (3) of the main text. The computation of ∆ pert a LO-HVP ,f (Q>Q max ) is described in Sec. 7. The second term in Eq. (5) is required to shift the renormalization point in Eq. (S5) from Q max to Q=0. It is obtained with lattice results for C f ii (t) through Eq. (4) of the main text, with Q=Q max . Its determination in the continuum limit is explained in Sec. 5.
The factor γ (Q max ) in Eq. (5) is simply Although calculating this function is trivial, for completeness we summarize its values for all three leptons and for the values of Q 2 max used here in Table S4.  In order for the separation of Eq. (5) to be valid, it must be independent of Q max within errors. We show that this is the case here. Since the perturbative contributions for a LO-HVP , =e, µ, are negligible, we focus on =τ . The dependence of a LO-HVP τ,f on Q max is plotted in Fig. S6 for f =ud, s, c, disc. As the figure indicates, the three terms of Eq. (5), for f =ud, disc, add up to a total a LO-HVP τ,f that is independent of Q 2 max within errors, for Q max ≥ √ 2 GeV. For f =s, c, where the lattice uncertainties are below one percent, this independence is no longer true within errors. However, the observed variations are accounted for by our matching uncertainty and are small compared to our total error on a LO-HVP τ .
The observed overall independence on Q max is an indication that our continuum-limit, lattice results are consistent with O(α 4 s ), five-loop perturbation theory for Q max ≥ √ 2 GeV. This result is highly nontrivial: it provides evidence that we control the high-virtuality contributions in our computation within our quoted errors. For Q < √ 2 GeV, agreement with perturbation is less good. Indeed, the results for Q max =1 GeV suggest that perturbation theory is beginning to break down at these low scales. tance quantities. For that purpose, we use the O(α 4 s ) result for the e + e − , R-ratio [57] and the optical theorem: where s is the center-of-mass energy, f in e + e − → ff X stands for a quark of flavor f and X can be gluons with, possibly, radiated quark pairs. For f =disc, which corresponds to the quark disconnected contribution, there are , which is the sum of the three contributions already mentioned, is represented by red squares. On these plots, error bars are smaller than the symbols, unless visible.
only gluon final states, which connect to the electromagnetic current through a quark loop. In the definition of R f pert , the denominator is the tree-level cross-section for e + e − → µ + µ − in the limit s m 2 µ . We then obtain the , needed in Eq. (S5), through the once subtracted dispersion relation: We compute this difference numerically, by feeding into the above dispersion relation, the four-loop output for R f pert (s) given by the routine rhad of [57]. The result for the difference is then redirected into Eq. (S5) and, again, integrated numerically. In rhad, we fix the strong coupling to its MS value at renormalization scale µ=Q max and use the physical values of the pole quark masses. Because the bottom and top contributions are negligible compared to our final errors and because we have not accounted for them in the lattice simulations from which we get the low Q 2 contributions, we neglect them here. The disconnected contributions only appear at O(α 3 s ) and vanish for the combined u, d and s terms, because the sum of these quarks' charges vanishes and because the effects of their masses are neglected in rhad. Moreover, the combined c, b and t quark-disconnected contributions are at least four orders of magnitude smaller than the perturbative corrections which we keep. Thus, we neglect disconnected contributions altogether. We do not include the tiny perturbative QED corrections, because they are also significantly smaller than our final errors and because we account for QED corrections separately as discussed below in Sec. 9. Finally, for Q max ≥ 1 GeV, the perturbative corrections to a LO-HVP e,f (Q≤Q max ) are all smaller than 10 −19 and can therefore be safely ignored.
The O(α 4 s ) results for ∆ pert a LO-HVP ,f (Q>Q max ) that we obtain are summarized in Table S5.
Only those for f =ud, c and =µ, τ are shown.
Because the u, d and s masses are neglected in rhad, which is reasonable given the Q max that we consider, we have ∆ pert a LO-HVP . All other perturbative corrections can be neglected.

Finite-volume corrections
As explained in the main text and in [58], the only finite-volume (FV) effects relevant in our percent-level (Q≤Qmax), =e, µ, and those for Π I=1 (Q 2 max ), which depend little on the momentum cut Qmax, in the range of interest, are given in the text.
calculation on lattices with spatial extents > ∼ 6 fm are assumed to be those associated with the two-pion contribution in the I=1 channel. We compute them to leading order in SU(2) chiral perturbation theory (χPT), as suggested in [58]. We study them for the scalar, A 1 representation of the cubic group that is appropriate for our determination of the HVP scalar function from the correlator 1 ii (t). We do not take into account discretization effects and, in particular, taste splittings, since we correct our continuum-extrapolated results.
We compute the corrections numerically. We first construct the lattice, position-space, pion propagator from the fast-Fourier-transformed, momentum-space, scalar propagator. Two such propagators are then appropriately combined to give the ππ contribution to the correlator C I=1 ii (t). The resulting correlator is subsequently treated in the same way as the corresponding lattice QCD, quark-antiquark correlator, to give the two-pion contribution to a LO-HVP ,I=1 (Q≤Q max ) andΠ I=1 (Q 2 max ). In particular, this means that the FV effects associated with the interpolation in Q, that is described in the paragraph preceding the "Lattice details" section of the main text, are accounted for. The bounding procedure is not implemented because uncertainties associated with this procedure are estimated independently.
To obtain these finite-volume corrections, we take L=6 fm and T =3L/2 which, within the uncertainties of these estimates, describes all of our lattices accurately. The infinite-volume result is taken to correspond to T = L = 13 fm and the relevant finite-volume differences are extrapolated to the continuum limit in the scalar theory to eliminate small discretization errors. Since we do not have the lattice simulations to check the accuracy of these predictions, we ascribe to them a 100% error and thus treat them as an order of magnitude estimate.
Because  − mu) corrections which must be added to our lattice results for the total LO-HVP contribution to the anomalous magnetic moments of the e, µ and τ to be able to compare them to those obtained from phenomenology. The individual corrections are described in the text. The total QED+(m d − mu) correction is given on the last line.
dependence. However, within the 100% uncertainty that we ascribe to these corrections, it is safe to take ∆ FV a LO-HVP µ,I=1 (Q≤Q max ) = 13.5 × 10 −10 for Q 2 max in the range 1 ÷ 5 GeV 2 . Only for those on a LO-HVP τ,I=1 (Q≤Q max ) is the Q max dependence significant. Those corrections are given in Table S6 for the Q max of interest. Note that these same, two-pion, FV corrections also affect the f =ud, disc contributions to a LO-HVP (Q≤Q max ). We have ∆ FV a LO-HVP (Q≤Q max ).
Our results forΠ I=1 (Q 2 max ) also suffer from FV corrections associated with the two-pion contribution. The dependence on Q 2 max of the corrections, that have to be added to our FV results forΠ I=1 (Q 2 max ), is mild in the range 1 ÷ 5 GeV 2 , with ∆ FVΠI=1 (Q 2 max ) = (2.3 ÷ 2.5) × 10 −4 . These same effects contribute toΠ ud (Q 2 max ) with a factor 10 9 and toΠ disc (Q 2 max ) with a coefficient − 1 9 . Here again we ascribe a 100% error to the estimates of these corrections.
Not surprisingly, the size of these FV corrections is larger, in relative terms, for the lighter leptons, since their anomalous moments are more sensitive to longerdistance physics. For the same reason, they are also slightly larger for smaller values of Q max , mostly in the case of the τ which is more sensitive to Q max .

QED and isospin-breaking corrections
Our calculation is performed in the isospin limit, with m u =m d and α=0. Moreover, as explained in the main text, we tune our lattice parameters to reproduce observables from an isospin-symmetric world in which m u =m d and α=0. Thus, our results for the LO-HVP contribution to the lepton anomalous magnetic moments are correct up to QED and (m d − m u ) effects that are proper to the HVP function. To compare our results with those determined using experimental cross sections, we must account for those isospin breaking (IB) effects. Eventually, this should be done ab initio with a lattice calculation. However, for the moment only exploratory calculations are available [37][38][39] and we resort here to phenomenology.
There are a variety of IB corrections to the lepton anomalous magnetic moments. The first, most easily quantifiable ones, are contributions from final states in the dispersive approach which are obviously absent in our calculation. These are the π 0 γ and ηγ contributions, which can be taken directly from the dispersive approach.
Here we take the results from [71].
A second set of IB contributions, which we are obviously missing, is final-state radiation (FSR). This can be determined by a combination of data and point-particle, QED corrections. Because of the positivity of the spectral function, this correction must be positive. Here we use the results from [71] and attribute to it a 50% error.
The third set of IB contributions requires hadronic models to estimate. The first of these is ρ-ω mixing, which we take from [71] and conservatively attribute to it a 50% error. The second has to do with other QED and (m d −m u ) effects on the hadrons which contribute to the HVP. One that is clearly identifiable is the fact that our charged pions have a mass of 134.8(3) MeV, which is lower than the physical M π ± =139.57 MeV. In the dispersive language, this means that in our calculation, the π + π − threshold is lower than it is in nature. The net effect is that our calculation overestimates the HVP contribution to the lepton anomalous magnetic moments. This effect is expected to be particularly pronounced in the case of the e and µ, which are sensitive to low energies. To estimate this effect, we take the difference between the π + π − contributions to a LO-HVP obtained to LO in χPT using M π ± and M π . We ascribe to it a 100% uncertainty to cover other, neglected effects.
In Table S7, we give quantitative estimates of these corrections to the LO-HVP contributions to the anomalous magnetic moments of all three leptons. The resulting total IB correction is given on the last line of Table S7 and included in our final results for the LO-HVP contributions to the anomalous magnetic moments of leptons of Table I from the main text.  Table I of the main text. In the results presented, the first error bar is statistical, the second is the systematic uncertainty associated with the continuum extrapolation, the third with the bounding procedure described in Sec. 2 (where applicable), the fourth with the matching to perturbation theory discussed in Sec. 6, the fifth with the lattice spacing uncertainty discussed in Sec. 3 and the sixth, where applicable, with FV corrections. Dashes in error brackets indicate that the corresponding systematic error does not affect the result in question.
from the main text, to obtain the individual flavor contributions to a LO-HVP , for =e, µ, τ . Thus, we take from Table S2 of Sec. 4 the lower vituality contributions, a LO-HVP ,f (Q<Q max ), obtained through the continuum limit of our lattice results. For the matching term, γ (Q max )Π f (Q 2 max ), we take the phase-space factor, γ (Q max ) of Eq. (S6), from Table S4. We obtain the nonperturbative quantity Π f (Q 2 max ) as described in Sec. 5 and take the results from Table S3. Finally, the high virtuality, perturbative contributions ∆ pert a LO-HVP ,f (Q>Q max ) are computed as explained in Sec. 7 and given in Table S5. To the f =ud, disc contributions, we have to add the FV corrections discussed in Sec. 8.
Our final results for the individual flavor contributions to the LO HVP component of the lepton anomalous magnetic moments are given in Table S8. These results include systematic errors associated with the continuum extrapolation, with our bounding procedure for the ud and disconnected contributions, with the matching to perturbation theory and with FV effects. These contributions are meant to be isospin limit quantities. Thus, we do not apply any QED or (m d − m u ) corrections to them. These are reserved for the total LO-HVP contri-bution, a LO-HVP , given in Table I of the main text.
In Fig. S7 we plot our results for the individual flavor contributions a LO-HVP µ,f , f =ud, s, c, disc, together with the only other ones available from the lattice [29,30,33,34,36]. Of those, only the results of [29,34,37] are obtained from N f =2 + 1 + 1 simulations. Those of [36] come from N f =2 and those of [30,33] from N f =2 + 1 simulations. For some reason, [30] do not include FV errors which are 2% of the I=0 contribution in our calculation, and should be at least as large in that reference.
As the figure shows, our ud contribution to a LO-HVP µ is significantly larger than the results of [34,36]. In particular, the difference with the only other N f =2 + 1 + 1 result published for this contribution is 2.2 combined standard deviations. Our result for the charm contribution is fully compatible with the two other lattice results. Finally, our result for a LO-HVP µ,disc is compatible with the only other determination [30] and, even with the inclusion of a FV uncertainty, it has a total error of 15%. This error represents 0.26% of a LO-HVP µ which means that it barely needs to be improved for determining a LO-HVP µ at the 0.2% level, as will be required by future experiments. with those from other lattice collaborations. The reference for the results plotted are HPQCD 14 [29], RBC/UKQCD 15 [30], RBC/UKQCD 16 [33], HPQCD 16 [34], Mainz 17 [36] and ETM [37]. For Mainz 17 latter, TMR and TMR+FV stand for the particular results of Table 7 in that paper.