Electromagnetic and strong isospin-breaking corrections to the muon $g - 2$ from Lattice QCD+QED

We present a lattice calculation of the leading-order electromagnetic and strong isospin-breaking corrections to the hadronic vacuum polarization (HVP) contribution to the anomalous magnetic moment of the muon. We employ the gauge configurations generated by the European Twisted Mass Collaboration (ETMC) with $N_f = 2+1+1$ dynamical quarks at three values of the lattice spacing ($a \simeq 0.062, 0.082, 0.089$ fm) with pion masses between $\simeq 210$ and $\simeq 450$ MeV. The results are obtained adopting the RM123 approach in the quenched-QED approximation, which neglects the charges of the sea quarks. Quark disconnected diagrams are not included. After the extrapolations to the physical pion mass and to the continuum and infinite-volume limits the contributions of the light, strange and charm quarks are respectively equal to $\delta a_\mu^{\rm HVP}(ud) = 7.2 ~ (2.5) \cdot 10^{-10}$, $\delta a_\mu^{\rm HVP}(s) = -0.0219 ~ (30) \cdot 10^{-10}$ and $\delta a_\mu^{\rm HVP}(c) = -0.0344 ~ (21) \cdot 10^{-10}$. At leading order in $\alpha_{em}$ and $(m_d - m_u) / \Lambda_{QCD}$ we obtain $\delta a_\mu^{\rm HVP}(udsc) = 7.1 ~ (2.9) \cdot 10^{-10}$, which is currently the most accurate determination of the isospin-breaking corrections to $a_\mu^{\rm HVP}$.

The muon anomalous magnetic moment a µ ≡ (g − 2)/2 is one of the most precisely determined quantities in particle physics. It is experimentally known with an accuracy of 0.54 ppm [1] (BNL E821) and the current precision of the Standard Model (SM) prediction is at the level of 0.4 ppm [2].
The discrepancy between the experimental value, a exp µ , and the SM prediction, a SM µ , corresponds to 3.7 standard deviations, namely a exp µ − a SM µ = 27.1 (7.3) · 10 −10 [3]. Since the above tension may be an exciting indication of new physics (NP) beyond the SM, an intense research program is currently underway in order to achieve a significant improvement of the uncertainties. The forthcoming g − 2 experiments at Fermilab (E989) [4] and J-PARC (E34) [5] aim at reducing the experimental uncertainty by a factor of four, down to 0.14 ppm, making the comparison of the experimental value a exp µ with the theoretical prediction a SM µ one of the most stringent tests of the SM in the quest of NP effects. On the theoretical side, the main uncertainty on a SM µ comes from hadronic contributions, related to the hadronic vacuum polarization (HVP) and light-by-light terms [6,7]. With the planned reduction of the experimental error, the uncertainty of the hadronic corrections will soon become the main limitation of this SM test.
The theoretical predictions for the hadronic contribution a HVP µ have been traditionally obtained from experimental data using dispersion relations for relating the HVP function to the experimental cross section data for e + e − annihilation into hadrons [8,9]. An alternative approach, proposed in Refs. [10][11][12], is to compute a HVP µ in Lattice QCD from the Euclidean correlation function of two electromagnetic (em) currents. In this respect an impressive progress in the lattice determinations of a HVP µ , which at leading order in α em is a quantity of order O(α 2 em ), has been achieved in the last few years [13][14][15][16][17][18][19][20][21][22][23][24].
With the increasing precision of the lattice calculations, it becomes necessary to include strong and em isospin-breaking (IB) corrections, which contribute to the HVP at order O(α 2 em (m d − m u )/Λ QCD ) and O(α 3 em ), respectively. In Ref. [20] a lattice calculation of the IB corrections to the HVP contributions due to strange and charm quarks, δa HVP µ (s) and δa HVP µ (c) 1 , was carried out using the RM123 approach [25,26], which is based on the expansion of the path integral in powers of the mass difference (m d − m u ) and of the em coupling α em . The quenched QED (qQED) approximation, which neglects the effects of sea-quark charges, was adopted and quark disconnected contractions were not included because of the large statistical fluctuations of the corresponding signals. The dominant source of uncertainty in the results of Ref. [20] was related to the em corrections to the renormalization constant (RC) of the local vector current, computed through the axial Ward-Takahashi identity derived in the QCD+QED theory.
In this work we present our determination of the IB corrections to the HVP contribution due to the light u-and d-quarks, δa HVP µ (ud), using the same methods and lattice setup adopted in Ref. [20] in the case of the strange and charm contributions. A preliminary result for δa HVP µ (ud) was presented in Ref. [27]. Thanks to a recent nonperturbative evaluation of QCD+QED effects on the RCs of bilinear operators performed in Refs. [28,29] we can update the determinations of δa HVP µ (s) and δa HVP µ (c) made in Ref. [20], obtaining a drastic improvement of the uncertainty by a factor of ≈ 4 and ≈ 6, respectively.
Within the qQED approximation and neglecting quark-disconnected diagrams the main results of the present study are: where the errors come from (statistics + fitting procedure), input parameters, chiral extrapolation, finite-volume and discretization effects. Thus, we confirm that the em corrections δa HVP µ (s) and δa HVP µ (c) turn out to be negligible with respect to the current uncertainties of the corresponding lowest-order terms a HVP µ (s) = 53.1 (2.5) · 10 −10 and a HVP µ (c) = 14.75 (0.56) · 10 −10 determined in Ref. [20]. In the case of the u-and d-quarks our finding (1) corresponds to about 1.2% of the lowest-order value a HVP µ (ud) = 619 (17.8) · 10 −10 obtained recently in Ref. [24].
The paper is organized as follows. In section II we describe the evaluation of the em and strong IB corrections to the light-quark HVP contribution at order O(α 2 em (m d − m u )/Λ QCD ) and O(α 3 em ) using the RM123 approach [25,26]. Details of the lattice simulations are collected in the Appendix A. In section III we describe the extrapolation to the physical pion mass and to the continuum and infinite volume limits. Finally, section IV contains our conclusions and outlooks for future developments.

II. ISOSPIN-BREAKING CORRECTIONS IN THE RM123 APPROACH
We adopt the time-momentum representation for the evaluation of the HVP contribution a HVP µ to the muon (g − 2), namely [30] a HVP where the kernel function K µ (t) is given by with m µ being the muon mass. In Eq. (5) the quantity V (t) is the vector current-current Euclidean correlator defined as where is the em current operator with q f being the electric charge of the quark with flavor f in units of the electron charge e, while ... means the average of the T -product over gluon and quark fields.
We will limit ourselves to the HVP contribution of the light u and d quarks, indicated by For sake of simplicity we drop the suffix (conn), but it is understood that in the following we refer always to quark-connected contractions only.
In the RM123 method of Refs. [25,26] the vector correlator V (t) is expanded into a lowestorder contribution V ud (t), evaluated in isospin symmetric QCD (i.e. m u = m d and α em = 0), and a correction δV ud (t) computed at leading order in the small parameters (m d − m u )/Λ QCD and α em : where the ellipses stand for higher order terms in (m d − m u )/Λ QCD and α em .
The separation between the isosymmetric QCD and the IB contributions, V ud (t) and δV ud (t), is prescription dependent. In this work we follow Ref. [20] and we impose the matching condition in which the renormalized coupling and quark masses in the full theory, α s and m f , and in isosym- s and m (0) f , coincide in the MS scheme at a scale of 2 GeV. Such a prescription is known as the GRS one [31].
The calculation of the IB correlator δV ud (t) requires the evaluation of the self-energy, exchange, tadpole, pseudoscalar and scalar insertion diagrams depicted in Fig. 1. More specifically one has where Z m is the mass RC in QCD only and the product Z m Z m encodes the corrections at first order in α em . The quantity Z m can be written as where Z (1) m is the pure QED contribution at leading order in α em , given in the MS scheme at the renormalization scale µ by [32,33] while Z f act m accounts for the corrections of order O(α n s ) with n ≥ 1 to Eq. (20). It represents the QCD correction to the "naive factorization" approximation Z m = Z ud being the renormalized light-quark mass in isosymmetric QCD.
In the numerical evaluation of the photon propagator the photon zero-mode has been removed according to the QED L prescription [34], i.e. the photon field A µ satisfies A µ (k 0 , k = 0) ≡ 0 for all k 0 .
In this work we make use of the same isosymmetric QCD gauge ensembles used in Ref. [20], i.e. those generated by the European Twisted Mass Collaboration (ETMC) with N f = 2 + 1 + 1 dynamical quarks, which include in the sea, besides two light mass-degenerate quarks, also the strange and the charm quarks with masses close to their physical values [35,36]. For earlier investigations of finite volume effects (FVEs) the ETMC produced three dedicated ensembles, A40.20, A40.24 and A40.32 (see Appendix A for details), which share the same light-quark mass and lattice spacing and differ only in the lattice size L. To improve such an investigation a further gauge ensemble, A40.40, has been generated at a larger value of the lattice size L.
For our maximally twisted-mass setup δm crit f has been determined in Ref. [37], while 1/Z m = Z P , where Z P is the RC of the pseudoscalar density evaluated in Ref. [38]. The coefficient Z f act m has been recently computed in Ref. [28] in a non-perturbative framework within the RI -MOM scheme [39].
Within the qQED approximation, which treats the dynamical quarks as electrically neutral particles, the correlator δV J (t) corresponds to the sum of the diagrams (1a)-(1b), while the correlators δV T (t) and δV P S (t) represent the contributions of the diagrams (1c) and (1d), respectively. The diagram (1e) contributes to both δV S (t) and δV SIB (t).
In our numerical simulations we have adopted the following local version of the vector current: where ψ f and ψ f represent two quarks with the same mass, charge and flavor, but regularized with opposite values of the Wilson r-parameter (i.e. r f = −r f ). Being at maximal twist the current (22) renormalizes multiplicatively with the RC Z A of the axial current. By construction the local current (22) does not generate quark-disconnected diagrams.
As discussed in Ref. [20], the properties of the kernel function K µ (t), given by Eq. (6), guarantee that the contact terms, generated in the HVP tensor by a local vector current, do not contribute to both a HV P µ and its IB correction.
Since we have adopted the renormalized vector current (22), the contribution δV Z A (t), appearing in Eq. (11), takes into account the em corrections to the RC Z A in QCD+QED, namely where Z A is the RC of the axial current in pure QCD (determined in Ref. [38]), while the product Z A Z A encodes the corrections at first order in α em . The quantity Z A can be written as where Z (1) A is the pure QED correction at leading order in α em , given by [32,33] Z (1) and Z f act A takes into account QCD corrections of order O(α n s ) with n ≥ 1 to Eq. (24). In this work we make use of the non-perturbative determination obtained in Ref. [28] within the RI -MOM scheme, which improves significantly the value Z f act A = 0.9 (1) obtained through the axial Ward-Takahashi identity in Ref. [20]. The values adopted for the coefficients Z f act m and Z f act A are collected in Table IV of Appendix A.
Thus, the IB term δV Z A (t) is simply given by where is the lowest-order contribution of the light-quarks to the vector correlator, calculated for our lattice setup in Ref. [24].
To sum up, the IB corrections δV ud (t) can be written as the sum of two (prescription dependent) contributions as where and δV SIB (t) is given by Eq. (16).
Within the qQED approximation, where the shift δm crit f is proportional to q 2 f [37], and neglecting quark-disconnected diagrams the QED correlator δV QED (t) is proportional to Using as inputs the experimental charged-and neutral-kaon masses the value m d − m u = 2.38 (18) MeV was determined in Ref. [37] at the physical point in the MS(2 GeV) scheme. Such a value is adopted in Eq. (16) for all gauge ensembles.
In Fig. 2 we show the dependence of both δV QED (t) and δV SIB (t) on the time distance t in the case of the ETMC gauge ensemble D20.48 (see Appendix A).

III. RESULTS
A convenient procedure [18,20,24] consists in splitting Eq. (5) into two contributions corresponding to 0 ≤ t ≤ T data and t > T data , respectively. In the first contribution the vector correlator is numerically evaluated on the lattice, while for the second contribution an analytic representation is required. If T data is large enough that the ground-state contribution is dominant for t > T data and smaller than T /2 in order to avoid backward signals, the IB corrections δa HVP µ (ud) can be written as with δa HVP where M ud V is the ground-state mass of the lowest-order correlator V ud (t) and Z ud V is the squared matrix element of the vector current between the ground-state |V and the vacuum: [24] the ground-state masses M ud V and the matrix elements Z ud V have been determined using appropriate time intervals t min ≤ t ≤ t max for each value of β and of the lattice volume for the ETMC ensembles adopted in this work. For the reader's convenience the values chosen for t min and t max in Ref. [24] are shown in Table I I. Values of t min and t max adopted in Ref. [24] to extract the ground-state signal from the light-quark vector correlator V ud (t) for each value of β and of the lattice volume V /a 4 for the ETMC gauge ensembles adopted in this work (see Table III of Appendix A).
In Eq. (31) the quantities δM ud V and δZ ud V can be extracted respectively from the "slope" and the "intercept" of the ratio δV ud (t)/V ud (t) at large time distances (see Refs. [20,25,26,37] where is almost a linear function of the Euclidean time t. This procedure is shown in Fig. 3 in the case of the gauge ensemble D20.48. The time dependencies of the integrand functions K µ (t) δV QED (t) and K µ (t) δV SIB (t) are shown in Fig. 4 in the case of the ETMC gauge ensemble B55.32 (see Appendix A). After summation over the time distance t, the SIB contribution dominates over the QED one.
The results for the separate contributions δa HVP  Table I), together with the uncertainty (at 1σ level) of a linear fit applied to the data.    (see Table I for the values of t min and t max , and it will not be given separately in the final error budget.
We have considered also the ratio of the IB correction δa HVP  For the separate QED and SIB contributions the FVEs differ qualitatively and quantitatively.
In the case of the QED data a power-law behavior in terms of the inverse lattice size 1/L is expected. According to the general findings of Ref. [40] the universal, structure-indepedent FVEs are expected to vanish, since they depend on the global charge of the meson states appearing in the spectral decomposition of the vector correlator, while the structure-dependent (SD) FVEs start at order O(1/L 2 ). Moreover, using the effective field theory approach of Ref. [41] one may argue that in the case of mesons with vanishing charge radius (as the ones appearing in the vector correlator) the SD FVEs start at order O(1/L 3 ) (see also Ref. [20]). In the case of the SIB correlator (16), since a fixed value m d −m u = 2.38 (18) MeV [37] is adopted for all gauge ensembles, an exponential dependence in terms of the quantity M π L is expected [42].
In Fig. 6 the data for the QED and SIB contributions to the ratio δa HVP Since the SIB data dominate over the QED ones, the FVEs for the ratio δa HVP For the combined extrapolations to the physical pion mass and to the continuum and infinite-2 We remind the reader that the lowest-order term a HVP µ (ud) has nonnegligible FVEs, which are exponentially suppressed in terms of MπL [42] (see Fig. 9 of Ref. [24]). volume limits we adopt the following fit ansatz: where the FVE term is estimated by using alternatively one of the fitting functions with B 0 and f 0 being the leading-order low-energy constants of ChPT and M 2 ≡ 2B 0 m ud . For the chiral extrapolation we consider either a quadratic (δ 1l = 0 and δ 2 = 0) or a logarithmic (δ 1l = 0 and δ 2 = 0) dependence. Half of the difference of the corresponding results extrapolated to the physical pion mass is used to estimate the systematic uncertainty due to the chiral extrapolation.
Discretization effects are estimated by including (D = 0) or excluding (D = 0) the term proportional to a 2 in Eq. (34). The free parameters to be determined by the fitting procedure are δ 0 , δ 1 , δ 1l (or δ 2 ), D and F (or F n ).
In our combined fit (34) the values of the free parameters are determined by a χ 2 -minimization procedure adopting an uncorrelated χ 2 . The uncertainties on the fitting parameters do not depend on the χ 2 -value, because they are obtained by using the bootstrap samplings of Ref. [38]. This guarantees that all the correlations among the lattice data points and among the fitting parameters are properly taken into account. The quality of our fitting procedure is illustrated in Fig. 7.
where the errors come in the order from (statistics + fitting procedure), input parameters of the eight branches of the quark mass analysis of Ref. [38], chiral extrapolation, finite-volume and discretization effects. In Eq. (36) the uncertainty in the square brackets corresponds to the sum in quadrature of the statistical and systematic errors.
Using the leading-order result a HVP µ (ud) = 619.0 (17.8) · 10 −10 from Ref. [24], we obtain our determination of the leading-order IB corrections to a HVP The above results show that the IB correction (37) is dominated by the strong SU (2)-breaking term, which corresponds roughly to ≈ 85% of δa HVP µ (ud).
Thanks to the recent nonperturbative evaluation of QCD+QED effects on the RCs of bilinear operators performed in Refs. [28,29] we can update the determinations of the strange δa HVP  [20]. Recently [23] in the case of the strange contribution the RBC/UKQCD collaboration has found the result δa HVP µ (s) = −0.0149 (32)·10 −10 , which deviates from our finding (40) by ≈ 1.6 standard deviations.
The sum of our three results (37), (40) and (41)  Recently, in Ref. [23] one QED disconnected diagram has been calculated in the case of the u-and d-quark contribution and found to be of the same order of the corresponding QED connected term.
Thus, we estimate that the uncertainty related to the qQED approximation and to the neglect of quark-disconnected diagrams is approximately equal to our QED contribution (38), obtaining δa HVP µ (udsc) = 7.1 (2.6) (1.2) qQED+disc · 10 −10 = 7.1 (2.9) · 10 −10 , which agrees within the errors with the recent determinations based on dispersive analyses of the experimental cross section data for e + e − annihilation into hadrons (see Ref. [3] and references therein).

IV. CONCLUSIONS
We have presented a lattice calculation of the isospin-breaking corrections to the HVP contribution of light quarks to the anomalous magnetic moment of the muon at order O[α 2 em (m d − m u )/Λ QCD ] in the light-quark mass difference and O(α 3 em ) in the em coupling. We have employed the gauge configurations generated by ETMC with N f = 2+1+1 dynamical quarks at three values of the lattice spacing (a 0.062 − 0.089 fm) with pion masses in the range M π 210 − 450 MeV and with strange and charm quark masses tuned at their physical values determined in Ref. [38].
The calculation of the IB corrections has been carried out adopting the RM123 approach of Refs. [25,26], which is based on the expansion of the lattice path-integral in powers of the small parameters (m d − m u )/Λ QCD and α em , which are both of the order of O(1%).
In this work we have taken into account only connected diagrams in which each quark flavor contributes separately. The leading-order em contributions to the renormalization constant of the local version of the lattice vector current, adopted in this work, have been evaluated using a recent nonperturbative calculation performed within the RI -MOM scheme in Refs. [28,29] . Thanks to that we have updated also the determinations of the strange δa HVP µ (s) and charm δa HVP µ (c) IB contributions made in Ref. [20], obtaining a drastic improvement of the uncertainties.
New QCD simulations with N f = 2+1+1 dynamical quarks close to the physical pion point [43] and the evaluation of quark-disconnected diagrams are in progress. and of the product M π L for the 16 ETMC gauge ensembles with N f = 2 + 1 + 1 dynamical quarks used in this contribution (see Ref. [38]) and for the gauge ensemble, A40.40 added to improve the investigation of FVEs. The bare twisted masses µ σ and µ δ describe the strange and charm sea doublet according to Ref. [51].
The central values and errors of the pion mass are evaluated using the bootstrap events of the eight branches of the analysis of Ref. [38]. The valence quarks in the pion are regularized with opposite values of the Wilson r-parameter in order to guarantee that discretization effects on the pion mass are of order O(a 2 µ ud Λ QCD ).
randomly. In the case of the light-quark contribution we have used 160 stochastic sources (diagonal in the spin variable and dense in the color one) per each gauge configuration.