Sz Borsanyi

Z Fodor

C Hoelbling

T Kawanai

S Krieg

Laurent Lellouch

R Malak

K Miura

K K Szabo

C Torrero

B C Toth

Hadronic vacuum polarization

come

Introduction.-Ever since the discovery of the electron's spin [START_REF] Gerlach | [END_REF]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 [START_REF] Bennett | Muon g-2)[END_REF] and of the standard model (SM) prediction (e.g. [START_REF] Davier | 14th International Workshop on Tau Lepton Physics (TAU 2016)[END_REF]) 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 [START_REF] Holzbauer | Proceedings, 12th International Conference on Beauty, Charm, and Hyperons in Hadronic Interactions (BEACH 2016)[END_REF], with the goal of reducing errors by a factor of four, and another one is planned at J-PARC, with similar objectives [START_REF] Otani | Proceedings, 2nd International Symposium on Science at J-PARC: Unlocking the Mysteries of Life, Matter and the Universe (J-PARC 2014[END_REF].

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 [START_REF] Fael | Proceedings, 12th International Workshop on Tau Lepton Physics (TAU 2012[END_REF] is being implemented [START_REF] Oberhof | Proceedings, 13th International Workshop on Tau Lepton Physics (TAU 2014)[END_REF].

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 78% of the total error [START_REF] Davier | 14th International Workshop on Tau Lepton Physics (TAU 2016)[END_REF]. This contribution also dominates the uncertainty in a τ [START_REF] Eidelman | [END_REF] and for a e , has a total size of around 6.5 times the experimental error [START_REF] Jegerlehner | Proceedings, Workshop on Flavour changing and conserving processes 2015 (FCCP2015[END_REF]. Today this contribution is determined most precisely using dispersion relations and the cross section of e + e -to hadrons and/or the rate of hadronic τ decays [3,[START_REF] Davier | 14th International Workshop on Tau Lepton Physics (TAU 2016)[END_REF][15][16][17][18][START_REF] Jegerlehner | KLOE-2 workshop on e + ecollider physics at 1 GeV, INFN-Laboratori Nazionali di[END_REF]. However, since the pioneering work of [START_REF] Blum | [END_REF], lattice QCD calculations of a LO-HVP µ [21][22][23][24][START_REF] Blum | Proceedings, 30th International Symposium on Lattice Field Theory (Lattice 2012)[END_REF][START_REF] Gregory | Proceedings, 31st International Symposium on Lattice Field Theory (Lattice 2013)[END_REF][START_REF] Malak | Proceedings, 32nd International Symposium on Lattice Field Theory (Lattice 2014)[END_REF][START_REF] Chakraborty | [END_REF][29][30][31][32][33][34][35][36][37][38] have made significant progress and provide a completely independent cross-check that will become competitive in the coming years.

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. [39] for a calculation of the b contribution to a LO-HVP µ ). Some previous lattice calculations of a LO-HVP , =e, µ [24,30], and of a LO-HVP , =e, µ, τ [40], 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 paper, we work directly at the physical point and compute disconnected contributions directly. Moreover, we implement a description of the lattice results [41,[START_REF] Spraggs | Proceedings, 33rd International Symposium on Lattice Field Theory (Lattice 2015)[END_REF] that solves the small virtuality issue [START_REF] Aubin | [END_REF] with finite-volume (FV) artefacts 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 arXiv:1711.04980v1 [hep-lat] 14 Nov 2017 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:

C µν (t) = 1 e 2 d 3 x j µ (x)j ν (0) , (1)
with e the positron charge, x=(t, x) and j µ /e= 2 3 ūγ µ u -1 3 dγ µ d -1 3 sγ µ s + 2 3 cγ µ 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,

C µν (t) = C ud µν (t) + C s µν (t) + C c µν (t) + C disc µν (t) = C I=1 µν (t) + C I=0 µν (t) , (2)
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 [START_REF] Blum | [END_REF]44]:

a LO-HVP ,f = α π 2 ∞ 0 dQ 2 m 2 ω Q 2 m 2 Πf (Q 2 ) , (3)
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 [41])

Πf (Q 2 ) ≡ Π f (Q 2 ) -Π f (0) (4) = 1 3 3 i=1 ∞ 0 dt t 2 - 4 Q 2 sin 2 Qt 2 Re C f ii (t) .
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 [START_REF] Malak | Proceedings, 32nd International Symposium on Lattice Field Theory (Lattice 2014)[END_REF] to be critical for reducing finite-volume (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:

a LO-HVP ,f = a LO-HVP ,f (Q≤Q max ) + γ (Q max ) Πf (Q 2 max ) +∆ pert a LO-HVP ,f (Q>Q max ) , (5)
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 [45]. 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 [45]. 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 [START_REF] Aubin | [END_REF], the Marichev interpolation advocated in [46] or the finite-energy sum rule approach of [47]. 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 [48] and a fermion action for N f =2 + 1 + 1 flavors of stout-smeared [49], 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 [50]. 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) [51] trajectories, are used. Details are given in [45] and more information about the simulations can be found in [52].

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 technique (AMA) of [53] 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 [29,54]. 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, as discussed in [34,45]. 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 (see [45]). 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. These fits are discussed in detail and examples are shown in [45]. Here we emphasize that with simulations at six lattice spacings down to 0.064 fm, we have full control over the continuum extrapolation.

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 [45], 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 ) [45]. 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 [55], as explained in [45]. 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 [45] we study the Q max dependence of our results for a LO-HVP . 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 three. 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 [45]. 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 [45].

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 [56,57], for large volumes those effects will be governed by pion contributions that can be computed in chiral perturbation theory (χPT) [56]. Since the I=0 channel is dominated by 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 [45] 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.4(13.4) × 10 -10 for the µ with negligible Q max dependence in the range of interest. For the τ they range from

= e [×10 -14 ] = µ [×10 -10 ] = τ [×10 -8 ] a LO-HVP ,I=1
156.9(2.4)(2.1)(0.0)(0.0) ( 3)(0.1)(0.0)(1.5)(4.6)(0.9) 711.0(7.5)(8.0)(0.2)(0.0)(5.5)(13.4)(5.1) 341.1(0.8)(2.0)(0.0)(1.0)(1.5)(1.9)(1.1) TABLE I. 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.

9.3(9.3)×10

-9 to 1.8(1.8)×10 -8 for Q max = 1÷ √ 5
GeV. We associate with these corrections a 100% uncertainty included in our error budget.

Compared to phenomenological determinations of a LO-HVP µ [START_REF] Davier | 14th International Workshop on Tau Lepton Physics (TAU 2016)[END_REF][15][16][17][START_REF] Jegerlehner | KLOE-2 workshop on e + ecollider physics at 1 GeV, INFN-Laboratori Nazionali di[END_REF], our m d =m u calculation without QED is missing isospin-breaking (IB) effects. These are detailed in [45]. Here we note that the corrections to be added are (1.7 ± 1.6) × 10 -14 for the e, (7.8 ± 5.1) × 10 -10 for the µ and (3.4 ± 1.1) × 10 -8 for the τ .

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.2%, 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 [24,33]. 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 [START_REF] Davier | 14th International Workshop on Tau Lepton Physics (TAU 2016)[END_REF]17,[START_REF] Jegerlehner | KLOE-2 workshop on e + ecollider physics at 1 GeV, INFN-Laboratori Nazionali di[END_REF].

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 µ [START_REF] Bennett | Muon g-2)[END_REF], assuming that all other SM contributions are unchanged, i.e. assuming no new physics (NP). Using the SM contributions summarized in [START_REF] Davier | 14th International Workshop on Tau Lepton Physics (TAU 2016)[END_REF], 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 phenomenological 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 [33] which is 1.9 σ away. A more detailed flavor-by-flavor comparison is given in [45]. However, our result is consistent with those from phenomenology within about one 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 con- with the only other two N f = 2+1+1 lattice QCD calculations [24,33] and with recent ones obtained from phenomenology [START_REF] Davier | 14th International Workshop on Tau Lepton Physics (TAU 2016)[END_REF]17,[START_REF] Jegerlehner | KLOE-2 workshop on e + ecollider physics at 1 GeV, INFN-Laboratori Nazionali di[END_REF]. 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. firm 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 [33,40] and only one concerning a LO-HVP τ [40]. The results in [40] are a LO-HVP e = 1.782(64)(86) × 10 -12 and a LO-HVP τ = 3.41(8)(6) × 10 -6 and in [33], a LO-HVP e = 1.779(39)×10 -12 . From the point of view of phenomenology, a dispersive analysis similar to the one implemented for the muon gives a LO-HVP e = 1.846(12) × 10 -12 [START_REF] Jegerlehner | Proceedings, Workshop on Flavour changing and conserving processes 2015 (FCCP2015[END_REF] and a LO-HVP τ = 3.38(4) × 10 -6 [START_REF] Eidelman | [END_REF]. Comparing these results to ours in Table I, we find the following. The result of [33] for a LO-HVP e displays a tension with ours which is slightly larger 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 [40] and for a LO-HVP τ , it is approximately three times more precise. In fact, our result for the latter is more precise than the phenomenological one.

M 2 π /f 2 π M 2(φ) π /f 2(φ) π 0.96 0.98 1.00 1.02 1.04 1.06 (M 2 K -0.5 * M 2 π )/f 2 π (M 2(φ) K -0.5 * M 2(φ) π )/f 2(φ) π β = 3.7000 (1) β = 3.7500 (3) β = 3.7753 (3) β = 3.8400 (3) β = 3.9200 (4) β = 4.0126 (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 β.

SUPPLEMENTAL MATERIAL 1. 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 [34]. 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

0 ≤ C ud ii (t) ≤ C ud ii (t c ) ϕ(t) ϕ(t c ) , (S1)
where

ϕ(t) = cosh [E 2π (T /2 -t)
] 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) 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 3 fm, for all three leptons. In our analyses, we consider the ranges t c = (3.000 ± 0.134) fm and t c = (2.866 ± 0.134) fm for the light connected, timelike correlators. For each t c in these intervals, we replace a LO-HVP ,ud (Q≤2 GeV), obtained directly from the lattice correlation function by the mean of its upper and lower bounds, and then average this mean over each t c interval. The two intervals are used to estimate the systematic uncertainty associated with our procedure, as explained in the main text. 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.

β a [fm] (T ×L/a 2 ) a 2 ∆taste #conf-ud/#conf-s/#conf-c/#conf-disc #src-ud/#src-s/#src-c/#src
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 ,

0 ≤ -C disc ii (t) ≤ 1 10 C ud ii (t c ) ϕ(t) ϕ(t c ) + C s ii (t) + C c ii (t) , (S2)
up to higher-order, wrap-around contributions. At large t, the connected strange and charm contributions in Eq. (S2) are exponentially suppressed, and their presence does not make a difference when determining t c , so we neglect them. As done in Fig. S2 for a LO-HVP ,ud (Q≤2 GeV), in Fig. S3 we plot the resulting upper and lower bounds on a LO-HVP ,disc (Q≤2 GeV) as a function of t c . Here, the t c time ranges are taken to be t c = (2.600 ± 0.134) fm and t c = (2.466 ± 0.134) fm, as suggested by the region of the merging of the bounds in the figure .  Pion-pion interactions change the smallest two-pion momentum from 2π/L in that channel. Using the model of [58] 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 ,ud and a LO-HVP ,disc by a small fraction of the statistical error.

Physical point and lattice spacing uncertainty

We define the physical mass point by using the isospin corrected pion and kaon masses, Mπ = 134.8(3) MeV and MK = 494.2(3) MeV, from [59], as well as the electromagnetically corrected η c mass, M ηc = 2.9863 [START_REF] Malak | Proceedings, 32nd International Symposium on Lattice Field Theory (Lattice 2014)[END_REF] GeV of [60]. 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 [61]. This yields a well defined physical point in the isospin limit. In intermediate steps of the analysis, we use the Wilson-flow-based [62] w 0 -scale [START_REF] Borsanyi | Budapest-Marseille-Wuppertal[END_REF].

As noted in [41,64], and generalized here for individual flavor contributions,

lim m →0 a LO-HVP ,f m 2 = 4π 2 3 α π 2 dΠ f (Q 2 ) dQ 2 Q 2 =0 . (S3)
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% 3) and ( 4), the correlator C ud ii (t), i = 1, 2, 3, is used up to and including t = tc in the sum over t in Eq. ( 4). For t>tc, C ud ii (t) is replaced by the upper and lower bounds of Eq. (S1). The red open circles correspond to the upper bound, obtained by assuming that the correlator falls off exponentially with the lowest two-pion energy, for t>tc. The grey closed circles are the lower bound obtained by setting C ud ii (t) to zero for t>tc. Points are slightly shifted around tc so that they do not overlap. The green band represents the average of the two bounds, with errors. The green crosses indicate the values of tc and the corresponding bound-averages that we consider for obtaining a LO-HVP ,ud

. All results are for an ensemble at β = 3.9200. error. Using f π = 130.50(1)(3)(13) MeV from the PDG [START_REF] Rosner | Submitted to: Particle Data Book[END_REF][START_REF] Patrignani | [END_REF], this means that we have a 0.4% error on our lattice spacings. 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 ,f (Q≤Q max ), are obtained by fits of our simulation results to

a LO-HVP ,f,lat (Q ≤ Q max ) = a LO-HVP ,f (Q≤Q max ) 1 +γ a,f a 2 + γ π,f (M lat π ) 2 -M 2 π +γ K,f (M lat K χ ) 2 -M 2 K χ (S4) +γ c,f M lat ηc -M ηc , with M 2 K χ = M 2 K -M 2 π /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 ,f,lat (Q≤2 GeV), = e, µ, τ and f =ud, s, c, disc, onto the a 2 -a LO-HVP ,f (Q≤2 GeV) planes, are shown in Fig. S4. The features of these fits are very similar for all three leptons, as they are for the other values of Q max (=1, √ 2, √ 3 and √ 5 GeV), so that we treat FIG. S4. Continuum extrapolation of the various flavor contributions to a LO-HVP (Q≤2 GeV) obtained using tc = (3.000 ± 0.134) fm for the ud contribution and tc = (2.600 ± 0.134) fm for the disconnected one. From left to right, = e, µ, τ . From top to bottom, the connected light, strange, charm, and disconnected contributions. The red open circles with errors are the results from our 15 simulations for the ud and s, 13 for the charm and 12 for the disconnected contributions, with statistical uncertainties. These points have been interpolated to the physical mass point using the fits to all lattice spacings (solid lines).

The different lines represent the fits obtained by imposing cuts in a (solid for no cut, dashed for a ≤ 0.118 fm, dotted for a ≤ 0.111 fm, dot-dashed for a ≤ 0.095 fm). The fact that a few of the lines do not appear to fit the red points is due to the dependence on other lattice parameters in those fits, which is slightly different from the one corresponding to the solid line. The green squares are the continuum extrapolated results for the given Qmax and tc, with statistical and continuum extrapolation errors only.

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.

(f, Q 2 max ) = e [×10 -14 ] = µ [×10 -10 ] = τ [×10 -8 ] (ud, 1 GeV 2 )
169.2(2.6)(2.3)(0.0) 627.1(7.5)(7.9)(0.1) 89.1(0.4)(0.8)(0.0) (s, 1 GeV 2 ) 13.6(0.0)(0.0)(-) 53.0(0.0)(0.1)(-) 9.2(0.0)(0.0)(-) (c, 1 GeV 2 ) 3.5(0.0)(0.0)(-) 14.4(0.0)(0. S2. Lattice results for the LO-HVP contribution to the anomalous magnetic moments of the e, µ and τ leptons obtained from Eqs. ( 3) and ( 4), where the integral over Q 2 in Eq. ( 3) is cut off at Q 2 max . This quantity is denoted a LO-HVP ,f

(Q≤Qmax) in Eq. ( 5), for the leptons = e, µ, τ and for the flavor contributions f =ud, s, c, disc, I=0, . Note that a LO-HVP ,I=1

= 9 10 a LO-HVP ,ud

. 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 "high-energy" 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.

Determination of

Πf (Q 2 max )
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). 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. TABLE S3. Values of the renormalized, scalar polarization function, Πf (Q 2 max ), as obtained on the lattice through Eq. ( 4), 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 [33]. 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.

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):

a LO-HVP ,f = a LO-HVP ,f (Q≤Q max ) +γ (Q max ) Πf (Q 2 max ) +∆ pert a LO-HVP ,f (Q>Q max ) ,
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.

∆ pert a LO-HVP ,f (Q>Q max ) ≡ (S5) α π 2 ∞ Q 2 max dQ 2 m 2 ω Q 2 m 2 Π f pert (Q 2 ) -Π f pert (Q 2 max ) ,
where the kinematic function, ω(r), is defined after Eq. ( 3). 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) with Q=Q max . Its determination in the continuum limit is explained in Sec. 5.

The factor γ (Q max ) in Eq. ( 5) is simply

γ (Q 2 max ) = α π 2 ∞ Q 2 max dQ 2 m 2 ω(Q 2 /m 2 ) . (S6)
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 2 max [GeV 2 ] γe(Q 2 max ) γµ(Q 2 max ) γτ (Q 2 max ) 1 7.
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.

Calculation of

∆ pert a LO-HVP ,f (Q>Qmax)
The contributions ∆ pert a LO-HVP ,f (Q>Q max ), defined in Eq. (S5), can be computed in QCD perturbation theory for Q max sufficiently large, because they are short distance quantities. For that purpose, we use the O(α The red open circles with errors are the results from our 15 simulations for the ud and s, 13 for the charm and 12 for the disconnected contributions, with statistical uncertainties. These points have been interpolated to the physical mass point using the fits to all lattice spacings (solid lines). The different lines represent the fits obtained by imposing cuts in a (solid for no cut, dashed for a ≤ 0.118 fm, dotted for a ≤ 0.111 fm, dot-dashed for a ≤ 0.095 fm). The fact that a few of the lines do not appear to fit the red points is due to the dependence on other lattice parameters in those fits, which is slightly different from the one corresponding to the solid line. The green squares are the continuum extrapolated results for the given tc, with statistical and continuum extrapolation errors only.

result for the e + e -, R-ratio [55] and the optical theorem:

ImΠ f pert (s) = - R f pert (s) 12π , R f pert (s) ≡ σ(e + e -→ f f X) 4πα 2 /(3s) , (
S7) where s is the center-of-mass energy, f in e + e -→ f f 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 only gluon final states, which connect to the electromag- on the choice of the value of Q=Qmax (from top to bottom, f =ud, s, c, disc). All three terms of Eq. ( 5), as well as their sum, are plotted as a function of Q 2 max for each f . The upside-down orange triangles correspond to a LO-HVP ,f (Q≤Qmax), as obtained on the lattice using Eqs. ( 3) and ( 4). The blue triangles correspond to the term γ (Qmax) Πf (Q 2 max ) of Eq. ( 5), with Πf (Q 2 max ) computed on the lattice using Eq. ( 4). The green circles represent the O(α 4 s ) perturbative contribution ∆ pert a LO-HVP ,f (Q>Qmax). For a LO-HVP τ,disc this contribution is not shown because it is negligible. a LO-HVP τ,f , 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. netic 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 difference

[Π f pert (Q 2 ) -Π f pert (Q 2 max )]
, needed in Eq. (S5), through the once subtracted dispersion relation:

Π f pert (Q 2 ) -Π f pert (Q 2 max ) (S8) = Q 2 -Q 2 max 12π 2 ∞ 0 ds 1 (s + Q 2 )(s + Q 2 max ) R f pert (s) .
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 [55]. The result for the difference is then redirected into Eq. (S5) and, again, integrated numerically. In rhad, we fix the strong cou- pling 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.

(f, Q 2 max ) = µ [×10 -10 ] = τ [×10 -8 ] (ud, 1, GeV
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 [56], the only finite-volume (FV) effects relevant in our percent-level calculation on lattices with spatial extents > ∼ 6 fm are (Q≤Qmax), that have to be added to our lattice results for a LO-HVP τ,I=1 (Q≤Qmax) to correct for the fact that our simulations are performed in a finite volume. The results for ∆ FV a LO-HVP ,I=1

Q 2 max [GeV 2 ] = τ [×10 -
(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.

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 [56]. 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 3 3 i=1 C I=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 by first constructing the ππ correlators, C I=1 ii (t), from the fast-Fourier transformed meson propagator on a very fine lattice where discretization errors are negligible. This correlator is then treated in exactly the same way as the corresponding quark 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 described in the paragraph preceding the "Lattice details" section of the main text, as well as with those pertaining to the upper and lower bound procedure of Sec. 2, are accounted for.

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. 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 of the very weak dependence of a LO-HVP e on Q max , the FV corrections to a LO-HVP e,I=1

(Q≤Q max ) are ∆ FV a LO-HVP e,I=1 (Q≤Q max ) = 4.6 × 10 -14 , independent of Q 2 max in the range 1 ÷ 5 GeV 2 . The FV corrections to a LO-HVP µ,I=1 (Q≤Q max ) have a very slight Q max 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.4 × 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 S7. QED and (m d -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.

= e [×10 -14 ] = µ [×10 -10 ] = τ [×10 -8 ] π 0 γ 1.
∆ FV a LO-HVP ,ud

(Q≤Q max )= 10 9 ∆ FV a LO-HVP ,I=1
(Q≤Q max ) and ∆ FV a LO-HVP ,disc

(Q≤Q max )= -1 9 ∆ FV a LO-HVP ,I=1
(Q≤Q max ).

Our results for ΠI=1 (Q 2 max ) also suffer from FV corrections associated with the two-pion contribution.

Here again the dependence of the FV corrections to ΠI=1 (Q 2 max ) on Q 2 max is very mild in the range 1÷5 GeV 2 and it is safe to take ∆ FV ΠI=1 (Q 2 max ) = 1.3 × 10 -4 , with 100% error, in that range. 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 . 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 [36][37][38] 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 [67].

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 [67] 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 [67] 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. We are now in a position to put together all of the ingredients of our calculation, according to Eq. ( 5), 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 S8. Final lattice results for individual flavor contributions to the LO-HVP componenents of the anomalous magnetic moments of the =e, µ and τ leptons. Results for a LO-HVP ,I=0,1 and a LO-HVP are given in 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. S5. To the f =ud, disc contributions, we have to add the FV corrections discussed in Sec. 8.

Sec. 7 and given in Table

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 contribution, 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 [START_REF] Chakraborty | [END_REF]29,32,33,35]. Of those, only the results of [START_REF] Chakraborty | [END_REF]33,36] are obtained from N f =2 + 1 + 1 simulations. Those of [35] come from N f =2 and those of [29,32] from N f =2 + 1 simulations. For some reason, [29] 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 [33,35]. 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 [29] 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 [START_REF] Chakraborty | [END_REF], RBC/UKQCD 15 [29], RBC/UKQCD 16 [32], HPQCD 16 [33], Mainz 17 [35] and ETM [36]. For Mainz 17 latter, TMR and TMR+FV stand for the particular results of Table 7 in that paper.

FIG. 1 .

 1 FIG. 1. Comparison of our result for a LO-HVP µ

  FIG. S2. Upper/lower bounds on a LO-HVP,ud (Q≤2 GeV) as a function of the time cut, tc. From top to bottom, the plots correspond to = e, µ, τ . In calculating a LO-HVP ,ud (Q≤2 GeV) from Eqs. (3) and (4), the correlator C ud ii (t), i = 1, 2, 3, is used up to and including t = tc in the sum over t in Eq. (4). For t>tc, C ud ii (t) is replaced by the upper and lower bounds of Eq. (S1). The red open circles correspond to the upper bound, obtained by assuming that the correlator falls off exponentially with the lowest two-pion energy, for t>tc. The grey closed circles are the lower bound obtained by setting C ud ii (t) to zero for t>tc. Points are slightly shifted around tc so that they do not overlap. The green band represents the average of the two bounds, with errors. The green crosses indicate the values of tc and the corresponding bound-averages that we consider for obtaining aLO-HVP

  FIG. S3. Same a Fig.S2, but for a LO-HVP ,disc (Q≤2 GeV) vs tc.

4 .

 4 Continuum limit and mass interpolation of a LO-HVP ,f,lat

2 ]

 2 FIG. S5. Continuum extrapolation of the various flavor contributions to Π(4 GeV 2 ) obtained using tc = (3.000±0.134) fm for the ud contribution and tc = (2.600 ± 0.134) fm for the disconnected one. From top to bottom, the connected light, strange, charm, and disconnected contributions. For the light and disconnected contributions, these correspond to the later of the two tc ranges. The red open circles with errors are the results from our 15 simulations for the ud and s, 13 for the charm and 12 for the disconnected contributions, with statistical uncertainties. These points have been interpolated to the physical mass point using the fits to all lattice spacings (solid lines). The different lines represent the fits obtained by imposing cuts in a (solid for no cut, dashed for a ≤ 0.118 fm, dotted for a ≤ 0.111 fm, dot-dashed for a ≤ 0.095 fm). The fact that a few of the lines do not appear to fit the red points is due to the dependence on other lattice parameters in those fits, which is slightly different from the one corresponding to the solid line. The green squares are the continuum extrapolated results for the given tc, with statistical and continuum extrapolation errors only.

  FIG. S6. Dependence of our results for a LO-HVP τ,f

10 ∆

 10 pert a LO-HVP ,ud + ∆ pert a LO-HVP ,c

10 .

 10 Individual flavor contributions to a LO-HVP and comparison with other lattice QCD calculations

  FIG. S7. Comparison of our results for the individual flavor contributions to a LO-HVP µ

TABLE S4 .

 S4 Values of the kinematical function γ (Qmax) defined in Eq. (S6), for =e, µ, τ .

  Q>Qmax for Q 2 max =1, • • • , 5 GeV 2 .Results are given for the connected ud and c contributions. All other contributions can either be obtained from those two or are completely negligible, as explained in the text.

	2 )	1.09	69.1
	(c, 1 GeV 2 )	0.12	11.5
	(ud, 2 GeV 2 )	0.26	34.9
	(c, 2 GeV 2 )	0.05	7.9
	(ud, 3 GeV 2 )	0.11	22.2
	(c, 3 GeV 2 )	0.03	6.1
	(ud, 4 GeV 2 )	0.06	15.6
	(c, 4 GeV 2 )	0.02	4.9
	(ud, 5 GeV 2 )	0.04	11.7
	(c, 5, GeV 2 )	0.01	4.1
	TABLE S5. Five-loop, O(α 4 s ) perturbative results for the contribution, ∆ pert a LO-HVP ,f (Q>Qmax) of Eq. (S5), to the
	anomalous magnetic moments of the leptons =µ and τ , com-
	ing from momenta

Acknowledgments.-We thank S. Dürr for contributions in the early stages of this project. In addition, LL thanks M. Benayoun, C. Davies, F. Jegerlehner, M. Knecht, C. Lehner, E. de Rafael and R. Van de Water for informative discussions. Computations were performed on JUQUEEN and JUROPA at Forschungszentrum Jülich, on Turing at the Institute for Development and Resources in Intensive Scientific Computing (IDRIS) in Orsay, on SuperMUC at Leibniz Supercomputing Centre in München, on Hermit at the High Performance Computing Center in Stuttgart. This project was supported by in part by the OCEVU Laboratoire d'excellence (ANR-11-LABX-0060) and the A*MIDEX Project (ANR-11-IDEX-0001-02), which are funded by the "Investissements d'avenir" French government program and managed by the "Agence nationale de la recherche" (ANR), by the DFG Grant SFB/TR55, by the Gauss Centre for Supercomputing e.V and by the GENCI-IDRIS supercomputing Grant No. 52275. RM was supported in part by a joint Ph.D. fellowship from the Centre national de la recherche scientifique (CNRS) and the Commissariat à l'énergie atomique et aux énergies alternatives (CEA).


Z Fodor

C Hoelbling

S Krieg

Laurent Lellouch

Th Lippert

A Portelli

A Sastre

K K Szabo

L Varnhorst

Up and down quark masses and corrections to Dashen's theorem from lattice QCD and quenched QED

published or not. The documents may come

The up (u) and down (d) quark masses are two fundamental parameters of the Standard Model of Particle Physics. These masses cannot be directly determined through experiment because of the confinement of quarks within hadrons. Lattice QCD provides an ab-initio approach to the non-perturbative calculation of QCD correlation functions. This method can be used to determine the light quark masses from the experimental values of hadron masses. In earlier work [2,3], we determined precisely m ud , the average of the up and down quark masses, using lattice QCD simulations at the physical values of the quark masses. This quantity has also been studied by many other lattice collaborations (cf. the FLAG review [4]) and considerable progress has been made on its determination. Thus, it is now relevant to aim for the calculation of the light-quark mass difference δm = m u -m d . This quantity is more difficult to obtain than m ud . Its small effect on hadron masses, of order O(δm/Λ QCD ) O(1%), is expected to be comparable in size to the leading O(α) electromagnetic (EM) corrections, usually not included in lattice simulations. Thus, earlier lattice calculations of this mass difference [2,3,5,[START_REF] Bazavov | Proceedings, 6th International Workshop on Chiral dymamics[END_REF] relied on phenomenological estimates of EM corrections. The inclusion of quenched QED effects was first performed in [START_REF] Duncan | [END_REF] on quenched QCD configurations, in [8,9] on N f = 2 QCD configurations and in [10] on N f = 2 + 1 configurations, at a single lattice spacing, with rather large pion masses and in small volumes. Preliminary results for the present calculation can be found in [11][START_REF] Lellouch | Light quarks and lattice QCD[END_REF][START_REF] Portelli | Proceedings, 32nd International Symposium on Lattice Field Theory (Lattice 2014)[END_REF] and preliminary results by MILC, in [START_REF] Basak | Proceedings, 26th Conference on Computational Physics (CCP 2014)[END_REF]. Very recently, an N f = 2 + 1 calculation of m u /m d in which QED effects are unquenched was presented in [START_REF] Horsley | [END_REF]. This calculation is performed on three ensembles at a single lattice spacing, in volumes up to (3.3 fm) 3 and with sea quark masses fixed at the SU(3) symmetric point M π = M K 412 MeV.

Here we include QED effects to the dynamics of the valence quarks, atop N f = 2 + 1 QCD configurations generated directly at the physical value of the light-quark masses, with full continuum and infinite-volume extrapolations, as well as with full non-perturbative renormalization and running. This is a sequel to the letter [1] which uses the same data set to compute light octet baryon isospin mass splittings. Note that a fully unquenched calculation of octet baryon and other hadron mass isospin splittings, with pion masses down to 195 MeV, can now be found in [START_REF] Borsanyi | [END_REF]. Here, because we are dealing with light quark masses whose extraction requires reaching deep into the chiral regime [17], we favor the simulations used in [1]. This data set also has the notable advantage that it has been used to determine the s and average u-d quark masses in [2,3]. Thus, all of the relevant non-perturbative renormalization and running has already been performed in pure QCD [2].

The light quark mass difference δm is connected, through a low energy theorem [18], to the pseudoscalar meson EM mass splittings. In the late 1960's, Dashen showed that pions and kaons receive the same EM contributions in the SU(3) chiral limit [19]. This result is commonly known as Dashen's theorem. During the 1990's, attempts to compute the chiral corrections to Dashen's theorem in effective field theories led to controversial and surprisingly large results (cf. the review [START_REF] Portelli | Proceedings, 32nd International Symposium on Lattice Field Theory (Lattice 2014)[END_REF] for more details). In this letter we present a computation of these corrections from our lattice QCD and quenched QED simulations.

General strategy. We consider in this work only the leading O(α, δm) corrections to isospin symmetry. As was done in [1], we define ∆M 2 to be the difference of the squared masses of the "connected" uu and dd pseudoscalar mesons. It is known from partiallyquenched chiral perturbation theory coupled to photons (PQχPT+QED) [20] that this quantity is related to δm by the following expansion:

∆M 2 = 2B 2 δm + O(m ud α, m ud δm, α 2 , αδm, δm 2 ) (1)
where B 2 is the two-flavor chiral condensate parameter. If the quark masses have their physical values, we can safely make the assumption that O(m ud ) = O(δm). Then at the level of precision considered here, ∆M 2 is proportional to δm. So to extract δm one needs to know the physical value of ∆M 2 and the constant B 2 .

B 2 was recently computed in [17], using the same QCD simulations as the ones considered in the present paper. To determine ∆M 2 , we consider the leading isospin expansion of the kaon mass splitting

∆M 2 K = M 2 K + -M 2 K0 : ∆M 2 K = C K α + D K ∆M 2 (2)
Results for ∆M 2 K obtained for different values of α and δm from lattice QCD and QED simulations can be fitted to this expression to obtain the coefficients C K and D K and subsequently the value of ∆M 2 corresponding to physical quark masses, from the experimental value of ∆M 2 K . Summary of the lattice methodology. The lattice setup used for this project is very similar to the one already described in [1]. The work is based on our set of lattice QCD simulations presented in [2]. It is composed of 47 N f = 2 + 1 QCD ensembles with pion masses down to 120 MeV, 5 lattice spacings down to 0.054 fm and 16 different volumes up to (6 fm) 3 . These simulations were performed using a tree-level O(a)-improved Wilson fermion action with 2 steps of HEX smearing. For each QCD configuration, a QED one is generated using the non-compact Maxwell action in Coulomb gauge with the four-momentum zero mode fixed to 0. The resulting SU (3) × U (1) configuration is then included in the Wilson-Dirac operator used to compute the valence quark propagators, with the appropriate electric charge. The valence light quark masses are tuned to explore the region where δm varies between 0 and its physical value. For most QCD ensembles, the unit of charge for valence quarks is set to its physical value. On one particular QCD ensemble, we perform three valence analyses: two 2). Here, ∆M 2 K is plotted as a function of ∆M 2 . The dependence of the lattice results on all other variables has been subtracted using the fit. The fit has a correlated χ 2 /dof equal to 1.59. It is plotted as a solid curve, with its 1σ band.

with close to physical δm and a value of α either about twice or one-fourth its physical value, and a third with α 0 and δm 0. A plot of the values of M 2 dd versus M 2 uu used in our valence datasets can be found in [1, Fig. 1].

In this setup, two approximations are made: the sea u and d quark masses have the same mass and they carry no electric charge (QED is quenched). It is straightforward to show that the splitting of the sea light-quark masses only affects isospin splittings at orders in the isospin expansion which are beyond those considered. Regarding the quenching of QED, large N c counting and SU(3) flavor symmetry suggests that the sea QED effects may represent O(10%) of the O(α) contribution to a given isospin splitting [1]. Considering the EM part of the kaon splitting, which is of particular interest here, the next-to-leading order (NLO) PQχPT+QED calculation of [20] can be used to estimate the QED quenching effects. In [21] we argue that they may represent 5% of the O(α) correction. Nevertheless, for giving the reader an idea of how such a quenching uncertainty may propagate to the other quantities studied in this paper, we retain the more conservative 10% quenching uncertainty on ∆ QED M 2 K = αC K .

The EM contribution to the kaon splitting. In the expansion (2), the coefficients C K and D K still depends on m ud , m s , a, and the temporal and spatial extents T and L. We fix m ud and m s to their physical values by matching to the experimental values of M 2 π + and the combination

(M 2 K + +M 2 K 0 -M 2 π + )/2.
Then, as explained in detail in [1], we use as a model for ∆M 2 K a first order expansion of C K and D K in these mass parameters around the physical mass point. Additionally, we allow for O(a) dis- cretization effects and power-like O(1/L) finite-volume effects [START_REF] Borsanyi | [END_REF]22] in the QED contribution proportional to C K , as required in our setup. The finite volume effects are taken into account by adding the following term to the aforementioned expansion of C K :

C FV K = - κM K L 1 + 2 M K L 1 - π 2κ T L + ρ L 3 (3)
where κ = 2.837 . . . is a known number, M K = 1 2 (M K + + M K 0 ) is the isospin averaged kaon mass and ρ is a fit parameter. In the parametrization (3), the O(1/L) and O(1/L 2 ) coefficients have been fixed to their known universal value [START_REF] Borsanyi | [END_REF] and the O(1/L 3 ) term is fitted to take into account additional structure-dependent effects. For the QCD contribution proportional to D K , we assume O(α s a, a 2 ) discretization effects and negligible finite-volume effects, which is justified in our large volumes given our present precision. To estimate systematic uncertainties, we consider a variety of analysis procedures. These variations are identical to those performed in [1]. They include (please see [1] for justifications and additional details): fitting the needed correlators on a conservative or a more aggressive time range; setting the scale with the mass of the Ω -or the isospin-averaged Ξ; eliminating points with M π + either greater than 400 MeV or than 450 MeV for the Ω -and the Ξ mass, and greater than 350 MeV or than 400 MeV for ∆M 2 K ; including either α s a or a 2 contributions in D K ; replacing individually the Taylor mass expansions in C K and D K by the inverse of these expansions (for a total of 4 choices). This leads to 128 different determinations of C K and D K . An example of such a fit is illustrated in Figs. 123. Finally, using the histogram method developed in [23], we combine all of these results to obtain:

∆ QED M 2 K = 2186(26)(68)(219) MeV 2 (4)
where C K is taken at the physical mass point, in the continuum and infinite volume limits. Here the first error is statistical, the second is systematic, and the third is an estimation of the quenching uncertainty as discussed above. Our result can be compared to an estimate obtained from the input of FLAG [4], ∆ QED M 2 K = 2090(380) MeV 2 . The results are entirely compatible and ours has a total precision which is more than 5 times higher omitting the generous estimate for the quenching error and more than 1.6 times including it. For completeness, we also give the value of the slope of ∆M 2 K in ∆M 2 at the physical point, obtained from our analysis: D K = 0.484(5)(4). This result is compatible with the value D K = 0.45 (9), obtained by appropriately combining results from FLAG [4]. Its total error is 15 times more precise.

Corrections to Dashen's theorem. As defined in [4], one can quantify corrections to Dashen's theorem with the parameter:

ε = ∆ QED M 2 K -∆ QED M 2 π ∆M 2 π (5)
The pion isospin mass splitting, needed to evaluate ε, is challenging to obtain through a lattice computation.

Because the neutral pion is diagonal in flavor, correlation functions for this state will contain quark disconnected diagrams. These diagrams are known to be expensive and hard to evaluate on the lattice. Thus, we choose not to compute the pion splitting here. Fortunately, using G-parity, one can easily show that the leading O(δm) corrections to ∆M 2 π vanishes. Therefore, at the level of precision considered in this paper, we have

∆ QED M 2 π = ∆M 2
π , which is very well known experimentally [24]. Using our result (4) for ∆ QED M 2 K and the experimen-tal value of ∆M 2 π , we obtain:

ε = 0.73(2)(5)(17) (6)
Here the relative quenching error, obtained by propagating a 10% uncertainty in ∆ QED M 2 K , is 23%. Now, if we include an estimate of the δm 2 corrections in the relation of ∆ QED M 2 π to ∆M 2 π , as given in [4] with the parameter ε m = ∆ QCD M 2 π /∆M 2 π = 0.04(2), we find ε = 0.77(2)(5) (17)(2), with the fourth uncertainty coming from the one in ε m . Our result of ( 6) can be compared to the FLAG estimate ε = 0.7(3) [4]. Again, it is entirely compatible with this estimate, and has a total precision which is more than 5 times higher without the quenching uncertainty estimate and about 1.8 times higher with it.

Up and down quark masses. Using our analysis of the kaon splitting, the experimental value of this splitting, our lattice result B 2 = 2.61 [START_REF] Bazavov | Proceedings, 6th International Workshop on Chiral dymamics[END_REF](1) GeV [17] in the MS scheme at 2 GeV and formula (2), we obtain:

δm = m u -m d = -2.41(6)(4)(9) MeV (7)
in the same scheme and at the same scale. If one assumes a 10% QED quenching error on ∆ QED M 2 K , this error propagates to 3.7% on δm. It is interesting to note that the quenching of QED has a rather small impact on the determination of This comes essentially from the fact that the QCD part of the kaon splitting is roughly 3 times larger than the QED Our result is entirely compatible with the value δm = -2.53 [START_REF] Borsanyi | [END_REF] MeV, derived from FLAG input [4]. Its precision is a factor of 2.1 to 1.3 higher, depending on whether the quenching error is taken into account.

If we combine [START_REF] Duncan | [END_REF] with our previous result m ud = 3.469(47)(48) MeV [3], we get: From the results of ( 8) and ( 9), we obtain the ratio of light quark masses:

m u = m ud + δm 2 = 2.
m u m d = 0.485(11)(8)(14) (10)
Q = m 2 s -m 2 ud m 2 d -m 2 u = 23.4(0.4)(0.3)(0.4) (11)
QED quenching effects of order 4% and 2%, respectively, cannot be excluded on these quantities. It is also interesting to compare our results for R to those obtained from χPT applied to η → 3π decays [18,[25][26][27][28]. The convergence of χPT for this process is very poor and it is usually supplemented by a dispersive analysis. Without such an analysis, the results vary from 19.1 at LO to 31.8 at NLO and 42.2 (or 38.7 setting the O(p 6 ) low-energy constants to 0) at NNLO [25]. The most recent NNLO χPT dispersive analysis [27] gives R = 37.7(2.2), in good agreement with our result. To summarize, our results are compatible with the estimates of [4], which already include input from the quenched QED studies mentioned above [29]. In most cases, they significantly improve on their precision. In all isospin symmetry breaking quantities the quenching uncertainty is the dominant one. Therefore, it is now important to determine these quantities using a fully unquenched calculation with significantly higher statistics, such as the one carried out in [START_REF] Borsanyi | [END_REF] for the splitting of stable hadrons. L.L. thanks Heiri Leutwyler and A.P. thanks Émilie Passemar for enlightening discussions. C.H., L.L. and A.P. acknowledge the warm welcome of the Benasque Center for Science where this paper was partly completed. Computations were performed using the PRACE Research Infrastructure resource JUGENE and JUROPA based in Germany at FZ Jülich, as well as HPC resources provided by the "Grand équipement national de calcul intensif" (GENCI) at the "Institut du développement et des ressources en informatique scientifique" (IDRIS) (GENCI-IDRIS Grant No. 52275), as well as clusters at Wuppertal and CPT. This work was supported in part by the OCEVU Labex (ANR-11-LABX-0060) and the A*MIDEX project (ANR-11-IDEX-0001-02), by CNRS grants GDR n 0 2921 and PICS n 0 4707, by EU grants FP7/2007-2013/ERC n 0 208740, MRTN-CT-2006-035482 (FLAVIAnet), by DFG grants FO 502/2, SFB TRR-55 and UK STFC Grants ST/L000296/1 and ST/L000458/1. L.V. was partially supported by a GSI grant.

Figure 1 .

 1 Figure 1. Example of a fit of the dependence of ∆M 2K on α and ∆M 2 to the expression of Eq. (2). Here, ∆M 2 K is plotted as a function of ∆M 2 . The dependence of the lattice results on all other variables has been subtracted using the fit. The fit has a correlated χ 2 /dof equal to 1.59. It is plotted as a solid curve, with its 1σ band.

Figure 2 .Figure 3 .

 23 Figure 2. Same fit as in Fig. 1. Here ∆QEDM 2 K is plotted as a function of 1/L. The greyed symbols show the full volume dependence of the data. For the plain symbols, the universal O(1/L) and O(1/L 2 ) finite volume effects from (3) have been subtracted and the fit to the O(1/L 3 ) correction is plotted.

  27(6)(5)(4) MeV (8) m d = m ud -δm 2 = 4.67(6)(5)(4) MeV (9) still in the MS scheme at 2 GeV. With the same assumptions as before, the QED quenching error on the individual quark masses is estimated to be 1.8 and 0.9%, respectively. Our results are nicely compatible with the FLAG values m u = 2.16(9)(7) MeV and m d = 4.68(14)(7) MeV.

  Our result is compatible with the FLAG average m u /m d = 0.46(2)(2). Moreover, our total precision is between 1.5 and 2.2 times higher, depending on whether our estimate of quenching uncertainties is included.We can further use our previous result m s /m ud = 27.53(20)(8)[3] to build the flavor breaking ratios R and

	Q:
	R =	m s -m ud m d -m u	= 38.2(1.1)(0.8)(1.4)

Strictly speaking, because u and d have different electric charges, this ratio is scale dependent in QCD plus QED. However it is easy to see that this dependency is beyond the leading isospin order considered in this work. Error propagations give a 2.9% QED quenching uncertainty on this ratio.