Strong-isospin-breaking correction to the muon anomalous magnetic moment from lattice QCD at the physical point

All lattice-QCD calculations of the hadronic-vacuum-polarization contribution to the muon's anomalous magnetic moment to-date have been performed with degenerate up- and down-quark masses. Here we calculate directly the strong-isospin-breaking correction to $a_\mu^{\rm HVP}$ for the first time with physical values of $m_u$ and $m_d$ and dynamical $u$, $d$, $s$, and $c$ quarks, thereby removing this important source of systematic uncertainty. We obtain a relative shift to be applied to lattice-QCD results obtained with degenerate light-quark masses of $\delta a_\mu^{{\rm HVP,} m_u \neq m_d}$= +1.5(7)%, in agreement with estimates from phenomenology and a recent lattice-QCD calculation with unphysically heavy pions.


I. INTRODUCTION
The anomalous magnetic moment of the muon, a µ ≡ (g µ − 2)/2, provides a stringent test of the Standard Model and a sensitive probe of new particles and forces beyond.It has been measured by BNL Experiment E821 to a precision of 0.54 parts-per-million [1], and the experimental result disagrees with Standard-Model theory expectations by more than three standard deviations [2].To investigate this discrepancy, the Muon g − 2 Experiment at Fermilab aims to reduce the experimental error by a factor of four, with a first result competitive with E821 expected within a year [3].To identify definitively whether any deviation observed between theory and experiment is due to new particles or forces, the theory error must be reduced to a commensurate level.
The largest source of theory uncertainty is from the O(α 2 EM ) hadronic vacuum-polarization contribution [2], a HVP µ , which is shown in Fig. 1.A well established method for determining this contribution employs dispersion relations combined with experimentally measured electron-positron inclusive scattering cross-section data.Recent determinations from this approach quote errors of 0.4-0.6%[4][5][6].Numerical lattice quantum chromodynamics (QCD) provides a complementary, systematic approach for calculating a HVP µ directly from first-principles QCD.Several independent lattice-QCD efforts to obtain a HVP µ are ongoing [7][8][9][10][11], with errors on recent determinations ranging from about 2-6% [8][9][10]12].The theoretical precision on a HVP µ needed to match the target experimental uncertainty is about 0.2%.
In this paper, we remove one of the largest systematic errors common to all current lattice-QCD calculations of a HVP µ , namely that arising from the use of degenerate up-and down-quark masses.We do so by calculating directly the strong-isospin-breaking contribution to the light-quark-connected contribution to a HVP µ .Phenomenological estimates suggest that the effect of strongisospin breaking on a HVP µ is about 1% [13][14][15].Electromagnetic effects are also not yet included in lattice-QCD calculations of a HVP µ and lead to a similar-sized uncertainty [16,17].In order to disentangle quarkmass from electromagnetic effects, we define the strong- isospin-breaking correction using up-and down-quark masses tuned from experimental hadron masses with QED effects removed [18].
The effect of strong-isospin breaking on the light-and strange-quark connected contributions to a HVP µ has been calculated in an exploratory study by the RBC/UKQCD Collaborations [19] in three-flavor lattice-QCD, with a heavy pion mass of about 340 MeV, and isospinsymmetric sea quarks.In this paper, we analyze two QCD gauge-field ensembles recently generated by the MILC Collaboration with four flavors of highly-improved staggered (HISQ) sea quarks and very high statistics; see Ref. [20] for methodology.One of the ensembles has fully nondegenerate quark masses with the u, d, s, and c quarks all fixed to their physical values.Our calculation is the first determination of δa HVP,mu =m d µ at the physical pion mass and with sea isospin breaking.
This Letter is organized as follows.We first present our numerical lattice-QCD calculation in Sec.II.Next, in Sec.III, we calculate the strong isospin-breaking correction to a HVP µ and discuss the contributions to the systemtatic error.We present our final result and compare it with phenomenological estimates and previous lattice-QCD calculations in Sec.IV.

II. LATTICE CALCULATION
We calculate the strong-isospin-breaking correction to a HVP µ on two new QCD gauge-field ensembles generated by the MILC Collaboration with four flavors of highlyimproved staggered (HISQ) quarks [20,21].Table I summarizes key parameters of the configurations.The two ensembles have the same approximate lattice spacing of a ≈ 0.15 fm, and the same strange-and charm-quark masses, which are both fixed close to their physical values.With staggered quarks, the pions possess an additional "taste" quantum number.Discretization errors from the HISQ action generate O(α 2 s a 2 ) corrections to the squared sea-pion masses of different tastes.On both ensembles, the mass of the taste-Goldstone ūd pion is fixed close to Nature's value of M π 0 ≈ 135 MeV, which is the mass that the charged pion would have in the absence of electromagnetism.The root-mean-squared pion TABLE I. Parameters of the QCD gauge-field ensembles.The first column shows the ratio of the lattice spacing to the gradient-flow scale w0 [22].To convert quantities in latticespacing units to GeV, we use w0 = 0.1715(9) fm [23].The next columns show the bare lattice up, down, strange, and charm sea-quark masses in lattice-spacing units and the number of configurations analyzed.The last column gives the taste-Goldstone sea-pion mass in GeV on each ensemble obtained from fits of pseudoscalar-current two-point correlators as in Ref. [20].Both ensembles have the same volume mass (averaging over tastes) is about 300 MeV.The two ensembles differ in one key feature: the values of the light sea-quark masses.Ensemble 1 is isospin symmetric, with the two light sea-quark masses equal and fixed to m l = (m u + m d )/2.Ensemble 2 features isospin breaking; here the two light-quark masses have the same average light-quark mass as ensemble 1, but the ratio of the light sea-quark masses is fixed to the value of m u /m d determined from the MILC Collaboration's study of pion and kaon electromagnetic mass splittings within the quenched approximation of QED [18].Comparing results on the two ensembles enables us to quantify the (tiny) effects of sea isospin breaking.
Our analysis strategy closely follows that of Ref. [8].On each ensemble, we calculate vector-current correlators j µ (x, t)j µ (0, 0) with zero spatial momentum and all four combinations of local and spatially-smeared interpolating operators at the source and sink.The smearing function is given in Eq. (A1) of Ref. [8], and we employ the same smearing parameters as in that work.To determine the quark-mass dependence of a HVP µ , we compute correlators with three valence-quark masses m q = {m u , (m u +m d )/2, m d }.With staggered quarks, the local vector current is not the conserved current, and must be renormalized.The renormalization factor Z qq V for HISQ quarks, however, has only mild quark-mass dependence so it cancels when the strong-isospin-breaking correction is calculated as a percentage shift.Additional details on the correlator construction and wavefunction smearings can be found in Ref. [8].
We fit the 2×2 matrix of correlators together using the multiexponential parameterization in Eq. (A6) of Ref. [8].The fit function includes contributions from both vector and opposite-parity states that arise with staggered valence quarks.The smeared correlators have smaller overlap with excited states than the local-local correlator, and therefore improve the determination of the energies and amplitudes.We fit the correlators over the symmetric time range [t min , T − t min ], thereby ensuring that the fit describes the correlator over the entire lattice time extent T .To reduce the degrees-of-freedom in the fit, in practice we average the correlator at times t and T − t and fit only to the lattice midpoint; we also average the smeared-source, local-sink correlator with the local-source, smeared-sink correlator.Because our limited number of configurations do not enable us to reliably determine the smallest eigenvalues of the correlation matrix, we employ singular-value-decomposition (SVD) cuts with the values chosen to obtain stable fits with good correlated χ 2 values.In practice, we replace all eigenvalues below the cut with the value of the SVD cut times the largest eigenvalue; this prescription increases the variance of the eigenmodes associated with the replaced eigenvalues and, thus, the errors on the fit parameters.We choose the number of states and fit range based on the stability of the ground-state and first-excited-state energies and amplitudes.
For both ensembles and all valence-quark masses, we obtain good correlated fits with stable central values and errors using t min /a ≥ 3, N states ≥ 3, and an SVD cut of 0.015, which modifies about 40% of the eigenvalues of the correlation matrix.For each of our six fits, the contribution to the χ 2 from the 66 correlator data points ranges from about 45 to 80.Although the lowest-energy states in the vector-current correlators are I = 1 ππ pairs, we do not see any evidence of such states in our two-point correlator fits.This is not surprising because there are only a few ππ states below the ρ mass in these correlators, and their amplitudes are suppressed by the reciprocal of the spatial volume.The ground-state energies for the correlators with m q = m l are E 0 = 776.7(6.5)MeV and E 0 = 779.4(5.1)MeV on the N f = 2 + 1 + 1 and N f = 1+1+1+1 ensembles, respectively; these are statistically consistent with the PDG average for the Breit-Wigner mass M ρ 0 = 775.26(25) MeV [24].
Following Ref. [8], we reduce the statistical errors in a HVP µ by replacing the correlator data at large times by the result of the multi-exponential fit.Although the fit function is appropriate for the periodic lattice temporal boundary conditions, we correct for the finite lattice temporal size by using the infinite-time fit function and doubling the correlator extent to t = 2T .We use the fitted correlator above t * > 1.5 fm; with this choice, roughly 80% of the value of a HVP µ comes from the data region.The values of a HVP µ computed with G fit (t) for t * > 1.5 fm agree within ∼ 1σ with those computed entirely from data, but with more than ten times smaller statistical errors for m q = m u .

III. ANALYSIS
We calculate a HVP µ using the method introduced by the HPQCD Collaboration [25], in which one constructs the [n, n] and [n, n − 1] Padé approximants for the renormalized hadronic vacuum polarization function ( Π(q 2 )) from time moments of zero-momentum vector-current correlation functions.These moments are proportional to the coefficients Π j in a Taylor expansion of Π(q 2 ) around q 2 = 0.The true result is guaranteed to lie between the [n, n] and [n, n − 1] Padé approximants.We employ the [3,3] Padé approximant for Π(q 2 ) obtained from the first six Taylor coefficients; the values of a HVP µ computed from the [3,2] and [3,3] Padé approximants differ by 0.1 × 10 −10 .
In Ref. [8], the [n, n] and [n, n − 1] Padé approximants for Π(q 2 ) are constructed from rescaled Taylor coefficients Π j × (E 0 /M ρ 0 ) 2j , where E 0 is the ground-state energy obtained from the two-point correlator fits.The rescaling was found to reduce the valence-quark-mass dependence of a HVP µ because the ρ-meson pole dominates the vacuum polarization.In addition, the rescaling cancels most of the error from the uncertainty on the lattice scale w 0 , which enters via the muon mass present in the one-loop QED integral for the a HVP,LO µ .Figure 2 shows a HVP µ on (1 + 1 + 1 + 1)-flavor ensemble at the up, down, and average light-quark masses.The valence-quark-mass dependence is statistically well resolved because the three points are strongly correlated, and is smaller after rescaling.
The physical value of the light-quark-connected contribution to a HVP µ is given by the sum of a HVP µ with two up quarks in the vector current and with two down quarks in the vector current weighted by the square of the quarks' electromagnetic charges: We then define the absolute shift with respect to the isospin-symmetric value as ∆a HVP,mu =m d µ ≡ a phys.
where m l ≡ (m u + m d )/2.Table II summarizes the isospin-breaking shifts on the N f = 2+1+1 and N f = 1+1+1+1 ensembles, both before and after rescaling the Taylor coefficients.The only significant source of systematic error in the shifts stems from the two-point correlator fits.The parametric errors from the lattice spacing are about a percent before rescaling, and are twenty times smaller with rescaling.The parametric errors from the current renormalization factor are ∼0.2%.The uncertainty due to the use of Padé approximants, which we take to be the difference between a HVP µ obtained with the [3,3] and [3,2] approximants, is about a percent.The 2.7% uncertainty on the ratio m u /m d in Ref. [18] stems largely from the estimate of electromagnetic effects, and leads to errors of about 2% and 1% on the physical up-and down-quark masses, respectively.Propagating the tuned quark-mass uncertainties to the physical a HVP µ using the measured slope of a HVP µ with respect to valence-quark mass changes the shifts in Table II by ∼ 0.2-0.3× 10 −10 ( < ∼ 0.1 × 10 −10 ) without (with) rescaling.Finally, the one-loop contributions from low-lying ππ states, which give rise to the leading finite-volume and discretization effects, cancel in ∆a HVP µ because the the charged pions in the loop are sensitive to the average of the up-and down-quark masses.
The shift in a HVP µ in Fig. 2 is smaller with rescaling because the valence-quark-mass dependence is milder.As expected, we do not observe any significant difference between the two ensembles.The leading sea isospin-breaking contributions to a HVP µ are quadratic in the difference (m d − m u ); taking Λ QCD = 300 MeV gives a rough power-counting estimate of their size as (m u − m d ) 2 /Λ 2 QCD ∼ 0.01%.We obtain our final results for the relative correction δa HVP,mu =m d µ by averaging the values on the two ensembles.The current renormalization factor drops out of the relative correction, and the errors from the lattice spacing, the up-and down-quark masses, and from higher-order finite-volume and ππ contributions are also further suppressed.

IV. RESULT AND OUTLOOK
We obtain for the relative strong isospin-breaking correction to the light-quark connected contribution to the muon g − 2 hadronic vacuum polarization δa HVP,mu =m d µ +1.5(4)% direct, (4.1) +0.4(4)% with E 0 rescaling, ( where the errors include Monte Carlo statistics and all systematics.Our result without rescaling the Taylor coefficients is compatible with the RBC/UKQCD Collaboration's calculation of the percentage shift using a pion mass of M π ≈ 340 MeV [19].It is also consistent with phenomenological estimates of the dominant isospinbreaking contribution from ρ-ω mixing using e + e − → π + π − data [13][14][15], ∆a ρ−ω mix.µ ∼ 2-5×10 −10 , and chiral perturbation theory [26], ∆a ρ−ω mix.µ ∼ 6 × 10 −10 , although ρ-ω mixing will also include effects from quarkline disconnected diagrams that we do not consider here. The percentage shifts in Eqs.(4.1) and (4.2) can be used to correct any existing or future result for the connected contribution to the hadronic vacuum polarization obtained with degenerate light quarks.Results for a HVP µ obtained without rescaling the Taylor coefficients should be corrected using Eq.(4.1); this applies to most recent lattice-QCD calculations.Equation (4.2) should be used to correct a HVP µ when E 0 rescaling is employed, as was done in the currently most precise lattice-QCD-calculation of a HVP µ from Ref. [8].Applying our result to 10 10 a HVP,ud µ = 599(6) stat.+sys.(8) QED+isospin from that work yields 10 10 a HVP,ud µ = 601(6) stat.+sys.(6) QED (2) isospin .With this reduction in the isospin-breaking error, the total uncertainty on the light-quark-connected contribution is limited by the ∼ 1% estimated error from the omission of the quarks' electromagnetic charges.We plan to directly calculate the electromagnetic correction to a HVP µ using dynamical QCD+QED gauge-field configurations to be generated soon by the MILC Collaboration with quarks, gluons, and photons in the sea [27].We anticipate that the inclusion of QED effects together with the isospinbreaking correction calculated here will enable a lattice-QCD calculation of the light-quark-connected contribution to a HVP µ at the 1% level.We have performed the first direct calculation of the strong-isosopin-breaking correction to a HVP µ at the physical up-and down-quark masses.We obtain an uncertainty on the relative correction of 0.4%, which is more than two times smaller, and also more reliable, than the ∼ 1% phenomenological estimate used in recent lattice-QCD calculations with equal up-and down-quark masses [8,10,28].Thus, it substantially reduces a significant source of uncertainty in a HVP µ , and is a crucial milestone towards a complete ab-initio lattice-QCD calculation of the hadronic contributions to a µ with the sub-percent precision needed by the Muon g − 2 and planned J-PARC experiments.

FIG. 1 .
FIG. 1.Leading hadronic contribution to aµ.The shaded circle denotes all corrections to the internal photon propagator from the vacuum polarization of u, d, s, c, and b quarks in the leading one-loop muon vertex diagram.

TABLE II .
Shift in a HVP µ from the isospin-symmetric to the physical valence-quark masses calculated on the ensembles in TableI.Results are shown both without and with rescaling the Taylor coefficients.As explained in the text, the numbers within a column should agree, but the two columns can (and should) differ.Errors shown include statistics and all systematic uncertainties.