Hadronic-vacuum-polarization contribution to the muon’s anomalous magnetic moment from four-flavor lattice QCD

SUPA, School of Physics and Astronomy, University of Glasgow, Glasgow G12 8QQ, United Kingdom Department of Physics and Astronomy, University of Utah, Salt Lake City, Utah 84112, USA Department of Physics, University of Illinois, Urbana, Illinois 61801, USA Fermi National Accelerator Laboratory, Batavia, Illinois 60510, USA CAFPE and Departamento de Física Teórica y del Cosmos, Universidad de Granada, 18071 Granada, Spain Department of Physics, Indiana University, Bloomington, Indiana 47405, USA Institute for Advanced Study, Technische Universität München, 85748 Garching, Germany Department of Physics, Syracuse University, Syracuse, New York 13244, USA Laboratory for Elementary-Particle Physics, Cornell University, Ithaca, New York 14853, USA Centre for Mathematical Sciences, University of Plymouth, Plymouth, PL4 8AA, United Kingdom Department of Physics, University of Colorado, Boulder, Colorado 80309, USA Department of Physics, University of Arizona, Tucson, Arizona 85721, USA


I. INTRODUCTION
In the absence of direct evidence for new particles or forces that are not present in the Standard Model, it becomes increasingly important to pursue experiments that may yield indirect evidence. Very heavy particles with masses beyond the reach of the Large Hadron Collider can have a tiny effect on low-energy observables through their brief appearance and disappearance in a quantum energy fluctuation of the vacuum that couples to the observable. Lighter particles with such small couplings to Standard-Model matter that they have escaped detection could behave in a similar way.
To pin down such effects requires both very precise experimental measurements and very good control of the theoretical calculations of the corresponding observables within the Standard-Model framework.
The anomalous magnetic moment of the muon, a μ , is such an observable. It is defined as ðg μ − 2Þ=2 where the gyromagnetic ratio, g μ , which connects the muon's spin and magnetic moment, would have a value of 2 in a world with no quantum corrections. Consequently, the value of a μ is sensitive to all of the particles that can appear virtually in a quantum-field-theory description of the muon/photon magnetic interaction. Given a careful enumeration of all of the Standard-Model contributions to a μ , we can identify any significant discrepancy with experiment as evidence for new physics.
The muon's anomalous magnetic moment was measured to an accuracy of 0.54 ppm nearly 20 years ago [1] at Brookhaven and will be updated to a planned accuracy of 0.14 ppm by the E989 experiment [2,3] now running at Fermilab and the E34 experiment [4] still under development at J-PARC. This prospect has galvanized a great deal of theoretical and linked experimental activity to improve the accuracy of the Standard-Model result for a μ . Recent calculations [5][6][7] give a Standard-Model result with an uncertainty at 0.3 ppm and a tantalizing 3.5-4 σ discrepancy with existing experiment. This theoretical precision is sufficient to achieve a greater than 5σ significance for the discrepancy if the central value does not change with the upcoming experimental results. It is nevertheless important to test the uncertainty in the Standard-Model result using different approaches to make sure that it is robust.
The results of Refs. [5][6][7] use experimental input for the cross section for e þ e − annihilation via a photon to hadrons as a function of center-of-mass energy to determine an important hadronic contribution to a μ known as the leading-order hadronic vacuum polarization contribution, a HVP;LO μ . This contribution, which appears at order α 2 , where α is the fine structure constant, is illustrated in Fig. 1. The uncertainty on its value is one of the two largest sources of error in the Standard-Model result. The leading-order hadronic vacuum polarization contribution can also be calculated from first principles using numerical lattice QCD, and there has been a great deal of progress in the past few years on improving lattice-QCD calculations of this quantity. 1 The aim of this effort is to reduce the uncertainty from lattice QCD first to a level commensurate with that from using σðe þ e − → hadronsÞ, and then to the ∼0.2% target precision of the Fermilab E989 and J-PARC experiments. In the meantime, however, lattice-QCD calculations already provide a strong test of those results from a completely different method with very different systematic errors.
As illustrated in Fig. 1, a HVP;LO μ requires knowledge of the quark vacuum-polarization function that couples to a photon [10,11]. In lattice QCD, individual diagrammatic contributions to the quark vacuum polarization can be considered separately via suitably constructed vector current-current correlation functions in Euclidean time. The vacuum polarization includes quark-line connected and disconnected diagrams, but the disconnected diagrams, where the quark loops are connected by intermediate gluons, contribute less than 2% to a HVP;LO μ [12][13][14][15][16][17][18]. The quark-connected contribution can be further separated into contributions from the individual quark flavors, up, down, strange, charm, and bottom. Accurate lattice-QCD results for the separate s-, c-, and (negligible) b-quark connected contributions to a HVP;LO μ were first obtained in Refs [19][20][21]. Subsequent lattice-QCD calculations [14][15][16][17][22][23][24] using different methods and quark formulations are in excellent agreement with these results.
The dominant quark-line connected contribution to a HVP;LO μ comes from the light (u=d) quarks, however, and is the target of this work. Here lattice-QCD calculations carry a number of additional technical challenges. The vector current-current correlator falls more slowly with Euclidean time at lighter quark masses, but at the same time the signal-to-noise degrades more rapidly. This means that the light-quark connected contribution to a HVP;LO μ receives contributions from larger Euclidean times than those from heavy quarks and that the data at these times are noisier. Hence, controlling statistical errors is a challenge. In addition, large physical volumes are needed for the lattice-QCD calculation to avoid systematic effects from squeezing light states (e.g., pions) into a small box.
The first lattice-QCD calculation of a ud μ ðconn:Þ, the light-quark connected contribution to a HVP;LO μ that included physical-mass u=d quarks was presented in Ref. [25], followed by several other lattice-QCD results [14][15][16][17]23]. FIG. 1. Leading hadronic contribution to the muon g μ − 2. 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. Diagrams in which the photon creates a quark-antiquark pair, which propagate while interacting via the strong and electromagnetic forces, and subsequently annihilate back into a photon, are called "quarkconnected" diagrams. Those in which the quark-antiquark pair annihilates into gluons are referred to as "quark-disconnected" diagrams.
All of these results were obtained in the isospin-symmetric limit, but the calculations differ in the quark formulation used, the lattice spacings and volumes available, and in the treatment of statistical errors and finite-volume effects. The agreement between different lattice-QCD calculations done independently will in the end be an important test of the results. Currently, the lattice-QCD results for a ud μ ðconn:Þ are spread over a range of several percent, with uncertainties at the same level. These errors are several times larger than those obtained using the experimental information from cross sections for e þ e − → hadrons. This means that lattice-QCD calculations are not yet in a position to add significant information to that available from e þ e − → hadrons [15]. This first round of complete lattice-QCD calculations has, however, crystallized the issues that must be addressed to improve current results and ultimately reach the target experimental precision.
In this paper, we present a calculation of the light-quark connected contribution to a HVP;LO μ in the isospin-symmetric limit. Like Ref. [25], our work uses the highly improved staggered quark (HISQ) action [26] and MILC ensembles with four flavors of HISQ sea quarks [27]. It also shares analysis strategies, a small set of common vector currentcurrent correlator data, and three coauthors with Ref. [25]. Many improvements have been made, however, with respect to that work. An important difference is that all ensembles of gluon-field configurations used in our analysis include u=d quarks of physical mass, eliminating the need for a chiral extrapolation, whereas in Ref. [25] only two out of ten ensembles were at physical u=d quark mass. In addition, we include ensembles at finer lattice spacings and one new ensemble that has approximately ten times the statistics of the others. The finer lattice spacings enable better control of the extrapolation to the continuum limit (zero lattice spacing), while the high-statistics ensemble allows us to undertake a significant study of the signal-tonoise issue mentioned above. This enables a better understanding of the impact of replacing correlator data with parametrizations of that data at large Euclidean times and will be discussed further in Sec. III A. There are also a number of differences in the analysis strategies employed in this work compared with Ref. [25], chiefly among them that the rescaling of the Taylor coefficients introduced in Ref. [25] is not used here. A detailed discussion of our analysis, including the differences with Ref. [25], is given in Secs. III B, III C, and IVA.
We do not present any new results for the contributions of strong-isospin-breaking and QED effects to the leadingorder hadronic vacuum polarization, nor for quark-line disconnected contributions. Progress has been made on all these small, but important, contributions recently [12][13][14][15]28,29]. We summarize the current situation for these pieces in Sec. IV to motivate the systematic uncertainty that we allow for not including them. This paper is organized as follows. Section II provides needed theoretical background to the calculation of the renormalized quark vacuum-polarization function from lattice QCD that is the key ingredient in calculating the leading-order hadronic vacuum polarization contribution to a μ . Section III gives details of our numerical lattice-QCD calculation and the methods we employ (along with datadriven tests of those methods) to tackle the issues of the growth of statistical uncertainties in the correlators and finite-volume effects. Section IV provides our results for the light-quark connected contribution to a HVP;LO μ and for the slope and curvature of the renormalized quark vacuumpolarization function, along with comprehensive error budgets for these quantities. Finally, Sec. V gives our determination of the total a HVP;LO μ from lattice QCD and compares it with other lattice-QCD results. This section also discusses the prospects for further improvements from lattice QCD that will allow significant input to be made to the Standard-Model value for a μ ahead of new experimental results.

II. BACKGROUND AND METHODOLOGY
The relation between the leading-order hadronicvacuum-polarization contribution to the muon's anomalous magnetic moment and the renormalized quark vacuumpolarization functionΠðQ 2 Þ ≡ ΠðQ 2 Þ − Πð0Þ, which is calculated here in lattice QCD, is given by [10,11] where Q denotes the Euclidean momentum carried by the virtual photons and K E ðQ 2 Þ is the standard kernel function introduced by Blum in Ref. [10]. The integrand peaks around Q 2 ≈ m 2 μ =2. The light-quark connected contribution to the muon's anomalous magnetic moment, a ud μ ðconn:Þ, arises from diagrams in which the photon in Fig. 1 produces light uū or dd pairs. We therefore start our lattice-QCD calculation of a ud μ ðconn:Þ with the zero-momentum u=d-quark currentcurrent correlation function in Euclidean space, where the summed index i runs over spatial components and j q i ¼qγ i q. The factor of 4 in front of the first term arises from the ratio of the quarks' electric charges squared, q 2 u =q 2 d . Following Ref. [19], we first compute time moments of GðtÞ, which are proportional to the coefficients Π j in a Taylor expansion ofΠðQ 2 Þ around Q 2 ¼ 0. We then obtain ΠðQ 2 Þ from ½n; n and ½n; n − 1 Padé approximants with n ¼ 3. BecauseΠðQ 2 Þ can be expressed in terms of a Stieltjes integral through a once-subtracted dispersion relation (see, e.g., Ref. [30]), the true result forΠðQ 2 Þ is guaranteed to lie between the ½n; n and ½n; n − 1 Padé approximants [31,32]. We find that the systematic uncertainty on a HVP;LO μ from the use of Padé approximants decreases with increasing n and is negligible even compared with the target experimental uncertainty for n ≥ 3. Indeed, we have checked with our lattice correlation functions that the time-moment method with n ¼ 3 Padé approximants as used in this work yields results for a ud μ ðconn:Þ that are numerically equivalent (to two decimal places or better) to the method introduced by Bernecker and Meyer in Ref. [33] based on the time-momentum representation of the Euclidean vector-current correlator.
In the time-momentum representation,Π is obtained directly from GðtÞ via the integral [33] which is a simpler procedure than calculating the Padé approximants from the time-moments of GðtÞ. However, the time-moment method directly yields the Taylor coefficients and hence allows us to correct them for finite-volume and lattice discretization effects using a chiral model of pions and ρ mesons before constructing a ud μ ðconn:Þ. In practice, then, the uncorrected values of a ud μ ðconn:Þ reported in Sec. III use Eq. (2.3) above, while the Taylor coefficients and corrected values of a ud μ ðconn:Þ are obtained from the time-moment method with n ¼ 3 Padé approximants.
The traditional and currently still most precise determinations of a HVP;LO μ use dispersive methods to obtain the vacuum-polarization function from experimental "R-ratio" data [5][6][7]34]Π where s is the square of the center-of-mass energy. With this approach, one integrates over all hadronic channels and it is not possible to cleanly identify which light-quark flavor was created at the photon vertex. Hence, one cannot separate their contributions to the cross section. One can, however, isolate heavy-quark contributions to the e þ e − cross section (see, e.g., Ref. [35]), enabling a clean comparison between lattice QCD and phenomenology. This is most clearly done at the level of the Taylor coefficients of the contribution toΠðQ 2 Þ for that quark flavor. The good agreement seen between lattice QCD c-and b-quark connected contributions to Π j and those from σðe þ e − → hadronsÞ [20,21,36] further substantiates the methods employed in the lattice-QCD calculations. In Sec. III C, we compare our lattice-QCD calculations of the Taylor coefficients Π j summed over all flavors with those from R-ratio data to check our model for calculating corrections due to nonzero lattice spacing and finite spatial volume.

III. LATTICE-QCD CALCULATION
We now present our lattice-QCD calculation. First, in Sec. III A, we describe the numerical simulations. We present the lattice quark and gluon actions employed and the parameters of the QCD gauge-field configurations and correlation functions. Next, in Sec. III B, we extract a ud μ ðconn:Þ in the isospin-symmetric limit on each ensemble from the vector-current correlation functions. We describe our approach for dealing with the substantial statistical noise in our two-point correlators at large times. Because we adapt many of the strategies of Ref. [25] in our analysis, we highlight key differences and improvements with respect to that work. Last, in Sec. III C, we correct the results for the isospin-symmetric a ud μ ðconn:Þ on each ensemble for finite-volume and taste-breaking discretization errors, and subsequently extrapolate these corrected values to zero lattice spacing.

A. Numerical simulations
We perform our calculation on QCD gauge-field configurations generated by the MILC Collaboration with four flavors of HISQ quarks [26,27]. These configurations are isospin-symmetric, i.e., the up and down sea-quark masses are equal with a mass m l ¼ ðm u þ m d Þ=2. We employ five ensembles with lattice spacings spanning a ≈ 0.15 − 0.06 fm and physical-mass light, strange, and charm sea quarks. The spatial volumes satisfy M π L ≳ 3.3 with M π the taste-Goldstone pion mass, while the temporal extents range (from coarsest to finest lattice spacing) between 7.2 ≳ T ≳ 10.2 fm. Table I summarizes key parameters of the configurations. Because our simulation light-quark masses are degenerate, throughout this work we use a ll μ ðconn:Þ to denote the quark-connected contribution from two light flavors in the isospin-symmetric limit. We reserve the notation a ud μ ðconn:Þ for nature's value. Two of the ensembles listed in Table I were also used in Ref. [25]: the a ≈ 0.15 fm ensemble with approximately 1000 configurations and the a ≈ 0.12 fm ensemble. Our analysis includes two new ensembles with a ≈ 0.09 fm and a ≈ 0.06 fm; the latter has a finer lattice spacing than those employed in Ref. [25], thereby providing better control over the continuum extrapolation. In addition, a new ensemble is included with a ≈ 0.15 fm and parameters identical to the older a ≈ 0.15 fm physical-mass ensemble, except for having better tuned quark masses. The new ensemble has 10 000 configurations, which is a factor of ten better statistics. On this ensemble, we can obtain a ll μ ðconn:Þ to high precision directly from the lattice vector-current correlator as described in Sec. II. Thus, comparing between this high-statistics ensemble and the older low-statistics one enables us to test our methods for extracting a ll μ ðconn:Þ from noisy data. Because we employ only physical-mass ensembles, a chiral extrapolation is not needed.
Following Ref. [25], on each ensemble we construct zero-momentum vector-current correlators with the valence-quark mass equal to the light sea-quark mass and four combinations of local and spatially smeared interpolating operators at the source and sink. We use the taste-vector current that combines quark and antiquark propagators at a single lattice site. The spatially smeared interpolating operators have the same taste because we employ a smearing function that combines separations of an even number of lattice spacings. This function is given in Eq. (A1) of Ref. [25], where the smearing parameters are also listed for lattice spacings a ≈ 0.15 − 0.09 fm. For the a ≈ 0.06 fm ensemble, we use a smearing radius that is the same in physical units as the one employed at a ≈ 0.09 fm, which yields the smearing parameters r 0 ¼ 6.75 and n ¼ 100. The correlators with smeared interpolating operators improve our identification of low-lying energy levels, to be discussed in Sec. III B. We take the correlators on the low-statistics a ≈ 0.15 fm and a ≈ 0.12 fm ensembles directly from Ref. [25]. These correlators were computed with 16 equally spaced random-wall time sources and averaged to gain statistics.
On the three newer ensembles analyzed in this work, we employ in addition a cost-effective variance-reduction technique called the truncated solver method (TSM) [40]. With this approach, on each configuration we compute a large number of "sloppy" correlators with a large relative error of 10 −5 at a small cost, and a single "fine" correlator with a small relative error of 10 −8 . We correct the average of the sloppy results using the difference between the approximate and precise solutions on a single source. In practice, we calculate sloppy propagators with all 48 time sources on the high-statistics a ≈ 0.15 fm ensemble, and from 16 time sources on the a ≈ 0.09 and 0.06 fm ones. Use of the TSM reduces our computational cost by more than a factor of 2.

B. Extraction of muon anomaly
A challenge common to all lattice-QCD calculations of a ll μ ðconn:Þ is the large statistical noise in the vector-current correlator at the physical light-quark mass, in particular for distances above about 2-3 fm. Figure 2 shows the locallocal vector-current correlator GðtÞ on the two a ≈ 0.15 fm ensembles. We average the correlator values at times t and T − t to increase statistics, and thus show the correlator only up to the lattice temporal midpoint. The low-and highstatistics data agree for times below 2 fm. Beyond this range, the data with low statistics become too noisy to yield a reliable estimate of the correlation function and hence of the contribution to a ll μ ðconn:Þ from large times. Several strategies to address the noise problem have been used in the literature [15,23,41]; here we follow the strategy of Ref. [25]. We first fit the 2 × 2 matrix of correlators with combinations of local and smeared sources and sinks together using the parametrization in Eq. (A2) of Ref. [25], constraining the energies and amplitudes with the Gaussian priors given in Eqs. (A3) and (A4) of Ref. [25]. In these fits, we minimize an augmented χ 2 that includes contributions from both the data and the priors [42]. Our fit function is simply a sum of exponentials expð−EtÞ such that the lowest-energy states are the only ones that survive to large time. With staggered quarks, the two-point correlators receive contributions from both correct parity and opposite-parity states; the latter lead to contributions that oscillate with time as ð−1Þ t . For every normal state in our fit (N states ), we also include an opposite parity state. We then replace the local-local correlator data for times above a chosen time t Ã by the result of the multiexponential fit, and TABLE I. Parameters of the QCD gauge-field ensembles. The first column shows the approximate lattice spacing, while the second lists the bare lattice up, down, strange, and charm sea-quark masses. The third column gives the ratio of the lattice spacing to the gradient-flow scale w 0 [37]; to convert quantities in lattice-spacing units to GeV, we use w 0 ¼ 0.1715ð9Þ fm [38]. The fourth column gives the nonperturbatively determined vector current renormalization factor obtained (for s quarks) in Ref. [39]. The fifth column lists the taste-Goldstone sea-pion masses; these were obtained from fits of pseudoscalar-current two-point correlators as in Ref. [27]. The sixth column shows the lowest-lying noninteracting two-pion energy level that couples to our vector current on each ensemble. The seventh column gives the lattice volumes. The final two columns give the number of configurations analyzed and the number of randomwall time sources used per configuration, where "TSM" indicates that we used the truncated solver method on this ensemble. use this mixed data þ fit correlator to calculate a ll μ ðconn:Þ either via Padé approximants or the time-momentum representation. Our detailed fit choices, e.g., fit ranges and number of states included, are given in Table II. They differ slightly from those of Ref. [25].
One must be careful with directly using the noisy largetime correlator data to calculate a ll μ ðconn:Þ. For all ensembles, we fix the maximum time (t max ) included in the fit based on plots of the local-local correlator (see Fig. 2), choosing t max slightly below the time at which the GðtÞ stops decaying exponentially. Beyond this point, the data violate the model-independent upper bound pointed out in Ref. [41] that GðtÞ must fall off more rapidly than expð−E ππ tÞ, where E ππ is the energy of two pions each with the smallest nonvanishing lattice momentum. The correlators stop decaying exponentially at around 2.3-2.6 fm on all ensembles with ∼1000 configurations. In contrast, the correlator on the ensemble with almost 10,000 configurations displays an exponential (cosh) fall-off until the lattice midpoint.
After fixing t max , we then vary the minimum time in the fit range (t min ) and the number of states in the fit function (N states ) and look for good correlated fits with stable central values and errors. Figure 3 plots a ll μ ðconn:Þ versus t min and N states on the a ≈ 0.09 fm ensemble. The inclusion of more states in the fit improves fits with smaller minimum times, and the a ll μ ðconn:Þ determinations are roughly independent of t min and N states for t min =a ≳ 8. The stability plots for other ensembles are qualitatively similar. Based on these plots, we choose t min ¼ 0.6 fm on the a ≈ 0.15 fm ensembles and increase t min smoothly with decreasing lattice spacing to t min ¼ 0.73 fm on the a ≈ 0.06 fm ensemble.
The true spectrum of the vector-current correlators is more complicated than the simple fit parametrization employed in our analysis, with many more levels than can be resolved within our finite statistics. Although we cannot identify the asymptotic lowest ππ energy level due to the large statistical noise in our data above around 2.5-3 fm, we can infer the presence of low-lying ππ states from the fitted ground-state energies, which are below the ρ pole on the finer ensembles. Even with these caveats, TABLE II. Parameters of the vector-current correlator fits. Ensembles are listed in the same order as in Table I. The number of degree-of-freedom is 3× the number of time slices in the fit range, rather than 4, because we average the local-source/ smeared-sink and smeared-source/local-sink correlators (which should be equal in the limit of infinite statistics) before fitting. The last column shows the standard frequentist p-values calculated from the χ 2 contribution from the data only (χ 2 data ) and with the degrees-of-freedom equal to the number of data points minus the number of fit parameters.  however, our fits provide a sufficiently accurate extrapolation of GðtÞ for the purposes of obtaining a ll μ ðconn:Þ. We have tested our noise-reduction strategy in several ways and summarized the studies that provide the strongest substantiation of our approach below.
The a ≈ 0.15 fm ensemble with ∼10 000 configurations enables a test of our use of correlator fits because on this ensemble we can obtain a ll μ ðconn:Þ reliably from data alone. Figure 4, left, shows the dependence of a ll μ ðconn:Þ computed from the mixed correlator on t Ã in fm for the two a ≈ 0.15 fm ensembles. Also shown is the 1σ error band for the value of a ll μ ðconn:Þ calculated entirely from data on the high-statistics ensemble. We find that, for all times t Ã ≲ 2.5 fm, the results for a ll μ ðconn:Þ obtained from G data ðt ≤ t Ã Þ and G fit ðt > t Ã Þ are consistent with the high-statistics data value. Further, the results on the lowand high-statistics ensembles are consistent with each other. This demonstrates that the fitted correlator yields an accurate value for a ll μ ðconn:Þ provided t Ã ≲ 2.5 fm. The number of low-lying ππ states in our vector-meson correlators increases rapidly as the lattice spacing, and consequently the taste splittings between sea-pion masses, decreases. Thus, it is also important to test our use of correlator fits with data that have several states below the ρ.
To obtain a correlator similar to our a ≈ 0.06 fm lattice data, but for which we know the spectrum exactly, we employ the chiral model in Appendix B of Ref. [25]. We first calculate the finite-volume energy levels, including ρ − ππ interactions, up to 2 GeV for our finest lattice spacing, a ≈ 0.06 fm. We then construct a fake correlator G fake ðtÞ with central values computed from the approximately 30 model energies and amplitudes, and a covariance matrix obtained from the simulation correlator G data ðtÞ. We then fit G fake ðtÞ using the same fit range as in our analysis and two or more states. Figure 5 plots G fake ðtÞ along with the result of a two-exponential fit. Figure 6 compares the individual contributions to a ll μ ðconn:Þ from each of the known states in G fake ðtÞ (top panel, blue) with those from each state in the two-state fit (bottom panel, red). Although the fitted energies are only a compromise between the actual energy levels, the value of a ll μ ðconn:Þ obtained from the fitted correlator (even with t Ã ¼ 0 fm) agrees with the known value to < 2 × 10 −10 . This is because the t-dependence of the fit correlator tracks G fake ðtÞ closely over the region of t that matters to a ll μ ðconn:Þ; the data are not sufficiently precise to distinguish between a two-state theory and the real theory. We have repeated this test using model spectra corresponding to each of our lattice spacings a ≈ 0.15 − 0.06 fm, and find the same conclusions. This indicates that our simple fit ansatz with two or more exponentials is sufficient to obtain the correct a ll μ ðconn:Þ to within the quoted statistics ⊕ fit uncertainties.
We also compare our approach with the bounding method used by the BMW Collaboration in Ref. [41]. With this approach, they select a value t c at which they replace the correlator data with the upper bound from a single exponential with the lowest-lying noninteracting two-pion energy level and a lower bound of zero. They then calculate a ll μ ðconn:Þ using the upper and lower bounds on the correlators varying the value of the matching point t c . They find that the upper and lower bounds meet at around 2.5-3 fm for their data, and take the average of a ll μ ðconn:Þ from the upper and lower bounds with t c ∼ 3 fm in their recent analysis [14]. Figure 7 compares a ll μ ðconn:Þ computed with our fit method and with BMW's bounding method on the low-statistics ≈0.15 fm ensemble. (See Table I for the relevant energy levels.) The results obtained with the two approaches agree, but the fit method yields smaller statistical errors on a ll μ ðconn:Þ. This is because a ll μ ðconn:Þ from the fit method is stable for t Ã above 1 fm, whereas the upper and lower bounds do not meet until around 2.5 fm, necessitating a larger value for t c . The consistency between the two noise-reduction strategies further substantiates our approach of using the fitted correlator at large times and also indicates that we obtain an accurate result for a ll μ ðconn:Þ with t Ã ∼ 1-2 fm. In Fig. 7, the value of a ll μ ðconn:Þ drifts upward beyond t Ã or t c around 2.5 fm, which corresponds to the time beyond which the correlator data no longer satisfy the modelindependent upper bound. We observe similar behavior on the other ensembles with only ∼1000 configurations. Thus, both the fit and bounding methods can overestimate a ll μ ðconn:Þ with noisy data if the replacement time t Ã or t c is chosen to be too large.
With the correlator fits in hand, we select the value of t Ã where we replace G data ðtÞ with G fit ðtÞ in our calculation of a ll μ ðconn:Þ. Plots of a ll μ ðconn:Þ versus t Ã show that the value of a ll μ ðconn:Þ is consistent within errors for t Ã between 0.5 and 2.5 fm. Our choice compromises between minimizing the statistical errors and maximizing the contributions from data. For simplicity, we select the same value of t Ã ¼ 2 fm for all ensembles, which is larger than the value t Ã ¼ 1.5 used in Ref. [25]. With our current choice, the data contribution to a ll μ ðconn:Þ is greater than 90% on all ensembles.

C. Lattice corrections and continuum extrapolation
Before we extrapolate the values obtained for a ll μ ðconn:Þ in Sec. III B to zero lattice spacing, we correct the data for the finite lattice spatial volume and for discretization effects from the mass splittings between staggered pions of different tastes. Both effects arise from one-loop diagrams with ππ intermediate states. As in Ref. [25], we calculate them within an extended chiral perturbation theory that includes pions, ρ mesons, and photons [43]. We work to one-pion-loop order, but to all orders in the leading interactions that couple the ρ 0 -γ-ππ channels. Details of the model calculation can be found in Appendix B in Ref. [25].
There are three differences between the numerical calculation of finite-volume corrections in Ref. [25] and in this work. The first difference is that the full one-loop finite-volume correction, which included a piece from quark-disconnected contributions, was applied to the raw a ll μ ðconn:Þ in Ref. [25]. Here we apply the quark-connected part of the one-loop finite-volume correction, which is 10=9 times the full one-loop value. Consequently our continuum-limit value of a ll μ ðconn:Þ will be larger than that in Ref. [25]. We address contributions to a HVP;LO μ from quark-disconnected contributions separately in Sec. IV B.
The second difference from Ref. [25] is that here we do not attempt to correct for differences between the simulated and physical values of the ρ meson's mass and decay constant by rescaling contributions to a μ . The majority of the lattice ensembles used in Ref. [25] have pion masses substantially larger than the physical pion mass; ρ rescaling was used to reduce the dependence of a ll μ ðconn:Þ on the light-quark mass. Here, however, all of our lattice ensembles use light-quark masses that are close to their physical values.
The third difference from Ref. [25] is a consequence of the second difference. Without ρ rescaling, we must include additional finite-volume corrections coming from the ρ's parameters (specifically Σð0Þ in Eqs. (B20) and (B22) in Ref. [25]). These new corrections are relatively small, adding 5 × 10 −10 to 13 × 10 −10 to a μ , depending upon the lattice ensemble. 2 For staggered quarks, the sea-pion masses are heavier than the taste-Goldstone pion for other representations of the approximate SO(4) taste symmetry. The taste splittings are discretization errors and thus decrease with lattice spacing. Consequently, the combined finite-volume plus discretization corrections are largest for our coarsest lattices, and decrease toward the continuum. The leading finite-volume correction to a ll μ ðconn:Þ in chiral perturbation theory is positive [44]. In total, the finite-volume plus discretization corrections to a ll μ ðconn:Þ for the lattice ensembles employed in our analysis range from approximately 68 × 10 −10 at a ≈ 0.15 fm to 31 × 10 −10 at a ≈ 0.06 fm. These include the leading-order contribution, from ππ loops, as well as next-to-leading-order corrections from the pion's charge radius and pion-pion scattering (see Appendix B of Ref. [25]). The next-to-leading-order corrections vary from ensemble to ensemble, but are are smaller than 6 × 10 −10 for our ensembles. Note that these subleading corrections are not included in the analyses of BMW [14] and RBC/UQKCD [15]. They are included, however, in the more recent analysis of Ref. [18]; our corrections are consistent with theirs (within errors).
As in Ref. [25], we can test our estimates of the lattice corrections by comparing our results for the Taylor coefficients of the vacuum-polarization function with phenomenological determinations from R-ratio data. Figure 8 compares our results for the total quark-connected contributions to Π 1 -Π 6 before and after the combined finitevolume plus discretization corrections are applied with a recent phenomenological determination by Keshavarzi et al. [7]. Because the experimental data include all possible diagrammatic contributions, for this test, we use the full one-loop correction, which includes both the connected and disconnected pieces. For our full range of lattice spacings, the corrections bring the lattice-QCD results into agreement with experiment, up to the 1%-2% level that might be expected from the small effects of strong-isospin breaking, QED, and quark-line disconnected diagrams missing from our calculation. Note that the high-n moments demonstrate that the continuum limit of our chiral theory agrees well with experiment, since the lattice contributions there are almost negligible (but these moments contribute little to a μ , as Fig. 8 also shows). These comparisons provide strong evidence that our estimated corrections are reliable both as a function of lattice volume and as a function of lattice spacing.
In Ref. [25], this model was also tested by comparison with an explicit finite-volume study on three a ≈ 0.12 fm ensembles with different spatial volumes but otherwise identical parameters. Because the pions were unphysically heavy on these lattices, there was little sensitivity to the spatial volumes. However, even the small spread in the raw results for a ll μ ðconn:Þ of 3(1)% was removed by the application of our combined finite-volume plus discretization corrections, providing further confidence in the method.
We can also compare our model with more recent finitevolume studies based on simulation results. These find finite-volume shifts of Δa ll μ ðconn:Þð5.4fm → 10.8fmÞ ¼ 40ð18Þ × 10 −10 , from the PACS Collaboration [45], and Δa ll μ ðconn:Þð4.66fm → 6.22fmÞ ¼ 21.6ð6.3Þ × 10 −10 , from the RBC/UKQCD Collaboration [46]. Our model, with all pion masses M π ¼ M π 0 and no staggered-pion mass splittings, gives shifts of 25ð4Þ × 10 −10 and 20ð3Þ × 10 −10 , respectively. These estimates agree with the lattice results above, within their large statistical uncertainties.  [7]. We take the s-, c-, and b-quark connected Π i s from HPQCD's companion calculations on the MILC HISQ ensembles [19,21]. 2 We only include the finite-volume part of this correction because effects due to the staggered-pion mass splittings have not been calculated (and vanish as a 2 → 0). Note also that we approximate parametersm andf by the physical ρ mass and decay constant, respectively, in the effective field theory used to calculate this correction (and all other finite-volume corrections); see Appendix B of Ref. [25]. Before extrapolating our results for a ll μ ðconn:Þ at nonzero lattice spacing to the continuum limit, we adjust the simulation values for the fact that our pion masses differ by a few MeV between ensembles (see Table I) and from the physical value. Using the same chiral model described above, we remove the continuum quark-connected contribution to a ll μ ðconn:Þ from γ → π þ π − → γ with the pion mass set equal to the simulation result for the Goldstone pion (and all other tastes of pion, once lattice artefacts are removed). We then reintroduce the continuum quarkconnected ππ contribution, but with the pion mass set equal to M π 0 ¼ 134.9766ð6Þ MeV [47]. Although the shifts are numerically tiny on the ensembles with M π 5 ∼ 135 MeV, the value of a ll μ ðconn:Þ on the outlying a ≈ 0.09 fm ensemble with M π 5 ∼ 128 MeV is decreased significantly, by about −8 × 10 −10 .
Finally, in order to account for higher-order contributions not included in the corrections, we assign 15% uncertainties to the net finite-volume and taste-breaking corrections on both a ll μ ðconn:Þ and the Taylor coefficients. Reference [25] assigned 10% uncertainties to these corrections. We use a larger uncertainty here because of the new sources of finite-volume error, associated with ρ parameters, that did not arise in the earlier analysis (see discussion above). These uncertainties are included in the errors on the corrected results listed in Table III and shown in Fig. 9. Figure 9 shows the lattice-spacing dependence of a ll μ ðconn:Þ before and after both lattice and M π corrections have been applied to the results obtained in Sec. III B, while Table III gives the numerical values. The net corrections range from about þ11% at a ≈ 0.15 fm to about þ5% at a ≈ 0.06 fm. Before corrections, the data display a large negative slope in a 2 . This is quite unlike what was seen for the s-quark connected contribution to a HVP;LO μ in Ref. [19], which also used the HISQ action and some of the same gauge-field ensembles as we use. There the variation with lattice spacing, from a ≈ 0.15 fm to the continuum, was only 0.5%. Most lattice-spacing artifacts are larger for squarks than for u=d-quarks, but taste-splittings are much larger for pions than for kaons. Hence the large latticespacing dependence seen here, before corrections are made, are almost certainly due to taste-splittings in the pion masses [25]. These should be greatly reduced by our corrections which account for the leading effects from taste splitting. Indeed, Fig. 9 shows no evidence at all of a 2 dependence in our corrected data. The fact that our combined finite-volume and discretization corrections remove the data's lattice-spacing dependence is perhaps the strongest evidence that our model for estimating these effects correctly describes the physics that underlies our numerical simulations.
We extrapolate the corrected values in Fig. 9 to the continuum limit using the following fit function, which TABLE III. Light-quark connected contribution to a HVP μ and the slope and curvature of the renormalized vacuum polarization before and after applying finite-volume, discretization, and M π corrections. Errors shown include uncertainties from statistics, two-point correlator fits, current renormalization, scale-setting, and finite-volume and discretization corrections. Results on all of the ensembles are correlated through the common scale w 0 used to convert m μ from physical to lattice units. Results on the two a ≈ 0.15 fm ensembles are also correlated to a much smaller extent because of the shared renormalization factor, Z V;ss .  where δm f ≡ m f − m phys f , and Λ ¼ 0.5 GeV is of order the QCD scale. This is similar to the fit function employed in Ref. [25], except that we no longer include terms to extrapolate in the valence-quark mass because all of our data are at the physical light-quark mass. The first term in parentheses adjusts for small sea-quark mass mistuning, while the second removes generic discretization errors; we employ priors for the coefficients c s ¼ 0.0ð3Þ and c a 2 ¼ 0ð1Þ. The values of a ll μ ðconn:Þ on each ensemble are statistically independent; we include in our fit correlations between the two a ≈ 0.15 fm ensembles from using the same Z V and between all ensembles from the common value of w 0 used to convert lattice-spacing units to GeV. Fitting our full data set to Eq. (3.1), we obtain a ll μ ðconn:Þ ¼ 637.8ð8.8Þ; c s ¼ 0.00ð30Þ; with a χ 2 =d:o:f: ¼ 0.04 and p ¼ 1. The fit posterior for both c s and c a 2 are consistent with zero, as expected because of the corrections applied to the data before extrapolation. Note that c a 2 ¼ −5ð1Þ for the raw values in Fig. 9.
To study the stability of the values and errors in Eq. (3.2), we consider a number of fit variations including adding higher-order terms in a 2 and δm f , doubling the prior widths on the fit parameters, and omitting the two coarsest ensembles. We show results for a ll μ ðconn:Þ for several of these variations in Fig. 10. Most variations differ only slightly from our original fit. The central values vary by no more than 16% of a standard deviation, while the uncertainties vary by at most 40% of a standard deviation. The fits are excellent, with χ 2 =d:o:f: < 0.1 in each case. The stability exhibited by these results suggest that our fit error accounts for the systematic uncertainties associated with the continuum extrapolation. The tiny χ 2 values suggest that our systematic errors are, if anything, overestimated.
We follow the same procedure to analyze the slope and curvature of the renormalized vacuum polarization, first applying finite-volume and taste-breaking discretization corrections, and then extrapolating to the continuum limit using Eq. (3.1). We obtain continuum-limit values of Π ll 1 ðconn:Þ ¼ 0.0932ð14Þ GeV 2 and Π ll 2 ðconn:Þ ¼ −0.2089ð64Þ GeV 4 . The p values of the fits are 1.0 and 0.8, respectively. The fit values for c a 2 are 0.02(87) and 0.23(98) for Π 1 and Π 2 , respectively; both are consistent with zero and similar to what we obtained for a ll μ ðconn:Þ. Also the sea-quark mass dependence of Π 1 and Π 2 is tiny, again like a ll μ ðconn:Þ. Finally, the continuum-limit values Π ll 1 ðconn:Þ and Π ll 2 ðconn:Þ are both stable against the fit variations discussed above for a ll μ ðconn:Þ.

IV. RESULTS
Here we present our final results for a ll μ ðconn:Þ, Π ll 1 , Π ll 2 , and a HVP;LO μ and the slope and curvature ofΠðQ 2 Þ with comprehensive error budgets.

A. Light-quark connected contribution
Our numerical calculation of a ll μ ðconn:Þ and the slope and curvature of the renormalized vacuum-polarization function described in the previous section is with equal up-and down-quark masses, and without electromagnetism. These corrections will be included a posteriori, as is done for other lattice-QCD g − 2 calculations in the literature. It is therefore useful to compare the available lattice-QCD results for a ll μ ðconn:Þ, Π ll 1 ðconn:Þ, and Π ll 2 ðconn:Þ before putting in the corrections for isospinbreaking and electromagnetism, in order to pin down the source of any disagreements among calculations. We employ the same definitions for the isospin limit of a ll μ ðconn:Þ, Π ll 1 ðconn:Þ, and Π ll 2 ðconn:Þ as in Refs. [14,15,23,48], which correspond to a world in which all pions have the same mass as the neutral pion. This allows for a clean comparison among lattice-QCD results. In Ref. [25], however, which appeared before Refs. [14,15,23,48], a different definition was used for a ll μ ðconn:Þ, which we describe below. Thus, the result of Ref. [25] for this quantity cannot be directly compared to ours or to those of Refs. [14,15,23,48].
Our results in the isospin-symmetric limit (taken from the fits in the previous section) are  We obtain a total uncertainty of 1.3% on the light-quark connected contribution to a HVP;LO μ in the isospin-symmetric limit without electromagnetism. The largest error contribution to Eq. (4.1) comes from the ∼0.5% uncertainty on the scale-setting parameter w 0 [38]. Because the Taylor coefficients of the vacuum-polarization function Π 1 has dimensions GeV −2 , the scale-setting error on a ll μ ðconn:Þ is approximately twice that of w 0 . Statistics, the continuum extrapolation, and finite-volume/discretization corrections also make significant contributions to the total error. The remaining contributions to the uncertainty in a ll μ ðconn:Þ are 0.1% or less.
In order to compare our result for a ll μ ðconn:Þ in Eq. (4.1) to the quantity reported in Ref. [25], we must account for the differences between definitions. Instead of quoting a value at the neutral pion mass as we do in this work, the a ll μ ðconn:Þ reported in Ref. [25] includes the one-loop continuum ππ contribution evaluated at the charged-pion mass. In addition, the corrections for finite volume and discretization effects applied in Ref. [25] include the quarkdisconnected contributions, while the corrections applied here include only the quark-connected contributions. The effects of both of these differences increase the value of a ll μ ðconn:Þ relative to Ref. [25]. After accounting for these differences, however, our result is still 1.8σ higher than the one in Ref. [25]. This is primarily because we do not rescale the Taylor coefficients by the ground-state energies of the correlator fits.
Despite the slightly different meanings of the light-quark connected contribution to a HVP;LO μ in Eq. (4.1) and in Ref. [25], it is still useful to compare the error budgets for these quantities. Compared with that work, we have reduced several key uncertainties. This is primarily because we employ only gauge-field configurations with physicalmass light quarks, two of which have finer lattice spacings than in that work. Consequently, the chiral extrapolation, which was an important source of error in Ref. [25], is replaced here by a chiral interpolation with an associated uncertainty of about 0.1%. Further, the error due to Padé approximants also made a significant contribution to the total uncertainty in Ref. [25]. It is reduced here to below 0.05% by using higher-order [3,2] and [3,3] Padés. Two of our uncertainty contributions in Table IV, however, are larger than in Ref. [25]. Because, in this analysis, we do not rescale the Taylor coefficients, our quoted lattice-spacing error is about 20 times larger than the estimate in that work. Our statistical and continuum-extrapolation errors are also two and three times larger, respectively, because the statistical errors increase with decreasing quark mass, and we only employ physical-mass light quarks. Overall, our total error on a ll μ ðconn:Þ is comparable to, but slightly larger than, the 1.1% error quoted in Ref. [25]. Note, however, that we have eliminated two systematic errors present in the result of Ref. [25] that were difficult to estimate and replaced them with statistical and systematic uncertainties that can be estimated more reliably. Figure 11 compares our result for a ll μ ðconn:Þ in Eq. (4.1) with recent unquenched lattice-QCD calculations [14-18, 23,48]. Our result is compatible with most of the IV. Error budgets for the Oðα 2 Þ light-quark-connected contribution, the leading Taylor coefficients of the vacuum-polarization function and the muon anomaly in the isospin-symmetric limit without electromagnetism. Sources of uncertainty that were considered, but found to have error contributions < 0.00%, are not shown. independently obtained values in the literature. Quantitatively, it agrees well with the published determinations by the BMW and ETM Collaborations [14,48], with the published results from Mainz (with N f ¼ 2) [23] and RBC/UKQCD [15], and with the calculation of Aubin et al. [18]. Our result for a ll μ ðconn:Þ is somewhat lower, however, than recent calculations (that appeared after this paper) by Mainz (with N f ¼ 3) [17] and Shintani and Kuramashi [16].
Finally we discuss the error budgets for the slope and curvature ofΠðQ 2 Þ, which are also given in Table IV. The uncertainty breakdown for Π ll 1 ðconn:Þ is similar to that for a ll μ ðconn:Þ because the two are proportional at lowest order in the Taylor expansion. The errors for Π ll 2 ðconn:Þ are different because it is more infrared than the other two quantities-the uncertainty due to uncalculated (higherorder) finite-volume/discretization contributions dominates all other contributions to the error budget. We do not quote values for higher-order Taylor coefficients ofΠðQ 2 Þ because the estimated errors from finite-volume plus taste-breaking discretization effects are no longer smaller than or commensurate with the contribution from statistics. Figure 12 compares our results for the slope and curvature of the renormalized vacuum-polarization function in Eqs. (4.2) and (4.3) with those from recent lattice-QCD calculations. Our result for the leading Taylor coefficient, Π ll 1 ðconn:Þ, agrees with those of the BMW [41] and ETM Collaborations [48], and is lower than that of the RBC/UKQCD Collaboration [15] by only 1.0σ. Our result for the second Taylor coefficient, Π ll 2 ðconn:Þ, agrees with the calculations of ETM and RBC/UKQCD, but is about 2.0σ larger in magnitude than that of BMW. The larger relative spread in Π ll 2 ðconn:Þ values between the collaborations may be due to the variety of approaches used to control the statistical error in the Euclidean vectorcurrent correlator at large times, since higher moments are sensitive to greater times.

B. Isospin-breaking, electromagnetic, and quark-disconnected contributions
To be able to compare our total summed over all quark flavors with experiment, we need to correct our result for a ll μ ðconn:Þ [Eq. (4.1)] for contributions due to strongisospin breaking, QED effects, and light-quark disconnected contributions. We will do this in four steps. First, we will consider these corrections for just diagrams with ππ intermediate states because they can be calculated reliably from the chiral model used in Sec. III C. Next, we will examine separately the remaining corrections from disconnected diagrams, strong isospin breaking, and QED. To estimate these contributions, we rely on our own lattice-QCD calculations when available, models, and phenomenology, and take generous uncertainties to cover roughly the spread of values in the literature. Table V [14][15][16][17][18]23,48]. All values correspond to isospin-symmetric QCD without electromagnetism. Results for a ll μ ðconn:Þ from four-, three, and two-flavor QCD simulations are denoted by squares, circles, and triangles, respectively. Note that the RBC/UKQCD Collaboration employed three-flavor QCD gauge-field configurations and then added the charm sea-quark contribution estimated from perturbation theory a posteriori.  [15,41,48]. All values correspond to isospin-symmetric QCD without electromagnetism. Note that we multiplied the Taylor coefficients quoted in Refs. [41,48] by the charge factor q 2 u þ q 2 d ¼ 5=9 so that they correspond to our normalization convention. estimated using the leading term in our chiral model. As discussed in Sec. III C, the chiral model gives an excellent description of the finite-volume and taste-breaking discretization effects in our numerical data and should therefore also be reliable here.
Because of spin-statistics, there is no π 0 π 0 contribution to a HVP;LO μ . Hence, the π 0 π 0 pieces must cancel between connected and disconnected diagrams. This leaves purely a π þ π − contribution, so it is clear that we should use the π þ mass when calculating corrections to our lattice-QCD result for a HVP;LO μ [25]. In Sec. III C, using our chiral model, we removed the continuum quark-connected contribution to a ll μ ðconn:Þ from γ → π þ π − → γ with the pion mass set equal to the simulation result for the Goldstone pion and then reintroduced it with the pion mass set equal to M π 0 . This is an artificial choice designed to yield a result for a ll μ ðconn:Þ in a world with equal u-and d-quark masses and without photons. Now, we can use our chiral model to subtract the continuum quark-connected ππ contribution with the pion mass set equal to M π 0 , and add the quark-connected contribution with the pion mass set equal to M π þ ¼ 139.57018ð35Þ MeV [44]. This yields for the part of isospin-breaking/electromagnetic correction coming from the π 0 π þ mass difference This correction already takes care of some QED effects because the difference between the π þ and π 0 masses comes largely from QED. Next, we calculate the contribution to a HVP;LO μ from quark-disconnected diagrams in Fig. 1 with ππ intermediate states. Because the ππ contribution appears only in the isospin-1 channel, the ratio of quark-disconnected to quarkconnected contributions is −1=10 from the ratio of appropriate quark electric charges [12,49]. Therefore, the ratio of quark-disconnected to total contributions is −1=9.
A calculation of the full ππ contribution to a HVP;LO μ within our chiral model using the experimental M π þ gives a ud μ ðππÞ ¼ 71 × 10 −10 [25]. Multiplying this by −1=9, we arrive at a quark-disconnected correction from ππ states of Δa ππ μ ðdisc:Þ ¼ −7.9 × 10 −10 : ð4:5Þ Adding Eqs. (4.4) and (4.5), we arrive at a total ππ correction to a HVP;LO μ from strong-isospin breaking, electromagnetism, and quark-disconnected diagrams of Δa ππ μ ¼ −12ð3Þ × 10 −10 : ð4:6Þ We assign a 25% error to this value because the dominant corrections to the leading-order ππ contribution in our chiral model (from the pion charge radius) enter at this level [25]. We follow the same prescription to estimate with our chiral model the ππ corrections to the slope and curvature ofΠðQ 2 Þ.

Residual light-quark disconnected corrections
There are also quark-line disconnected corrections to a HVP;LO μ that have nothing to do with the ππ contribution discussed above. Following the approach of Chakraborty et al., we estimate these by examining the contributions to the anomaly from the ρ and ω mesons [12]. Together, these two resonances account for almost 80% of the total a HVP;LO μ [5][6][7]. The ratio of the disconnected to connected moments coming from the ρ and ω is given by Eq. (11) in Ref. [12], where the moments (now) include the quarks' electric charge factors. This relation, when combined with experimental data for ρ and ω masses and bounds on their widths, implies a disconnected contribution from non-ππ states of Δa ρω μ ðdisc:Þ ¼ −5ð5Þ × 10 −10 ; ð4:8Þ where the error is from the uncertainty on the inputs. The correction in Eq. (4.8) does not include disconnected diagrams that mix light-quark and s-quark loops (connected to the photons), but these are known to be much smaller [12]. Again, we estimate the disconnected contribution from the ρ and ω resonances to the Taylor coefficients Π 1 and Π 2 in the same manner. Note that adding the above −5ð5Þ × 10 −10 to the ππ contribution from Eq. (4.5) gives −13ð5Þ × 10 −10 for the total quark-line disconnected contribution. This is well in line with direct lattice-QCD calculations of the quark-disconnected contribution to a HVP;LO μ in the isospin-symmetric limit and without QED-including s-quark contributions, the BMW  [15]-and further supports the reliability of our model calculations.

Residual strong-isospin breaking corrections
The effects from QCD-isospin breaking (i.e., quark-mass differences) and QED are intertwined both in nature and in lattice-QCD simulations because QED contributions shift the bare quark masses. Here we define the residual strongisospin correction to a HVP;LO μ as the shift relative to the isospin-symmetric value a ll μ ðconn:Þ that results when the bare u and d quark masses are retuned separately so that (i) their average gives the experimental value for the π 0 mass [as required for a ll μ ðconn:Þ], and (ii) their ratio has the physical value obtained from lattice-QCD calculations including electromagnetism [50,51]. Note that ππ contributions largely cancel in this correction because the pion mass is primarily sensitive to the average light-quark mass.
When only the quark-line connected diagrams are considered, the strong-isospin breaking correction will contain unphysical effects from ππ states where the π meson is composed of uū and dd states. These effects will be positive since isospin-breaking effects are positive and the "π u " meson is unnaturally light. They will be canceled, as discussed above, when the quark-line disconnected diagram is included. This means that we might expect substantial negative contributions from the quarkline disconnected diagrams, relative to the isospin-symmetric case, when strong-isospin breaking effects are included. Indeed, our preliminary results for the strongisospin-breaking correction to the quark-disconnected contribution confirm this [53]. We therefore increase the errors on our initial estimate of the total residual correction from strong-isospin breaking (from [28]) to allow for disconnected contributions of a commensurate size, giving Δa ud μ ðSIBÞ ¼ 10ð10Þ × 10 −10 : ð4:9Þ The analysis in Ref. [28] also yielded estimates for the strong-isospin breaking corrections to the Taylor coefficients ofΠðQ 2 Þ of δΠ HVP;LO 1 ðm u ≠ m d Þ ¼ þ1.6ð6Þ% and δΠ HVP;LO 2 ðm u ≠ m d Þ ¼ þ3.0ð8Þ%. We employ these values to obtain the absolute corrections to Π HVP;LO 1 and Π HVP;LO 2 and again increase the uncertainties to 100% to allow for large quark-disconnected contributions.

Residual QED corrections
We have already included a sizeable part of the full QED correction by replacing the π 0 mass by the π þ mass in the ππ contribution. We estimate the residual corrections from QED, beyond those accounted for above, via powercounting to be of order α ∼ 1%. This yields an estimate for the absolute correction to a HVP;LO μ of Δa ud μ ðQEDÞ ¼ 0ð5Þ × 10 −10 ; ð4:10Þ where we have taken a central value of zero because we do not know the sign of the correction. We take the same relative QED error for the Taylor coefficients Π 1 and Π 2 . Our estimate of residual QED corrections is consistent with results from the analysis of a HVP;LO μ based upon experimental data on e þ e − → hadrons. For example, the contribution from the simplest photon channel, e þ e − → π 0 γ, is 4.5 × 10 −10 [7]. Equation (4.10) is also consistent with (still early) efforts to estimate the QED contribution using lattice-QCD simulations [15,29,52]. The RBC/UKQCD Collaboration finds Δa HVP;LO μ ðQEDÞ ≈ −1ð6Þ × 10 −10 from summing results from connected and disconnected diagrams [15], while the ETM Collaboration finds Δa HVP;LO μ ðQED;conn:Þ ¼ 1.3ð1.0Þ × 10 −10 from connected diagrams only [52].

Total contribution from u=d quarks
Summing the corrections from Eqs. (4.6) and (4.8)-(4.10), we obtain for the total correction from strong-isospin breaking, QED, and quark-disconnected contributions where the first error is from a ll μ ðconn:Þ and the second is from Δa ud μ .

C. Total leading-order contribution
Finally, to obtain the total leading-order hadronic vacuum polarization contribution to a μ , we add the contributions from heavy flavors to a ud μ (Eq. (4.12). We take the connected results for strange, charm, and bottom quarks calculated by the HPQCD Collaboration Ref. [19][20][21]. 3 Disconnected contributions from these quarks are expected to be negligible compared with our other uncertainties. We follow the same procedure for the Taylor coefficients of the renormalized vacuum-polarization function. . More than 90% of the central value comes from the light-quark connected contribution, as does about 30% of the error. The remainder of the error on a HVP;LO μ comes from the uncertainty on our estimate of the missing contributions from QED, strongisospin breaking, and quark-disconnected diagrams. The contributions from s, c, and b quarks generate the remaining ∼10% of the central value, while contributing a negligible amount, ∼0.1%, to the error.

V. SUMMARY AND OUTLOOK
Our results for the leading-order HVP contributions to a μ and the slope and curvature of the renormalized vacuum polarization function are (Table VI)  The total uncertainty on a HVP;LO μ is ∼2.2% and is dominated by our conservative estimate of the combined uncertainty from the omission of strong isospin-breaking, electromagnetism, and quark-disconnected contributions in our calculation of the u=d-quark contribution (see Sec. IV B).
We also reiterate the key intermediate result of this work, which is our new determination of the light-quark connected contribution to a HVP;LO μ in the isospin-symmetric limit and without electromagnetism [from Eq. This result improves upon and supersedes the calculation of Chakraborty et al. in Ref. [25]. As can be seen from Fig. 11, our determination of a ll μ ðconn:Þ has smaller errors than other recent unquenched lattice-QCD calculations [14][15][16][17][18]23,48]. This is primarily because our fit method for controlling the statistical errors in the Euclidean vectorcurrent correlator at large times yields smaller uncertainties on a ll μ ðconn:Þ than approaches used by other collaborations. In addition, we include higher-order contributions in our calculation of the finite-volume corrections. (Several tests of this our fit method are summarized in Sec. III B and of our model for finite-volume corrections in Sec. III C.) Figure 13 compares our determination of the total, leading-order hadronic-vacuum-polarization contribution to a μ in Eq. (5.1) with other lattice-QCD calculations [14][15][16][17]23,48,54] and phenomenological analyses of experimental R-ratio data [5][6][7]55]. Our result agrees with all but one of the independent lattice calculations and has a comparable error. 4 It also agrees with the R-ratio analyses, although with roughly five to seven times larger uncertainties.
We also compare our result for a HVP;LO μ in Eq. (5.1) to the expectation from experiment. Assuming that there are no contributions to the muon anomalous magnetic moment from physics beyond the Standard Model, the BNL E821 Experiment [1] implies a value for a HVP;LO μ of 720ð7Þ × 10 −10 . This value is obtained by subtracting from experiment the calculated values of QED [56], electroweak TABLE VI. Individual flavor contributions to the leading Taylor coefficients of the vacuum-polarization function and the muon anomaly. The first error quoted for the u=d contributions is from the lattice analysis; the second comes from uncertainties in our estimates of the effects of strong isospin-breaking, electromagnetism, and quark disconnected diagrams. Results for strange and heavier quarks include only the quark-connected contributions and are not new, but come from earlier HPQCD calculations [19][20][21]; disconnected contributions are expected to be negligible. The definitions of the Taylor coefficients include the factor of the quark's electric charge squared.  3 The present author list overlaps with those of Refs. [19][20][21]. 4 As we were finishing this paper, Shintani and Kuramashi presented a new determination of 10 10 a HVP;LO μ ¼ 737ðþ16; −21Þ [16] that is 1.5σ above our result and is in more than 2σ-tension with the R-ratio analyses. [57], and higher order HVP [58,59] contributions and the consensus value for the hadronic light-by-light term [60]. Our result is 1.3σ below the "no new physics" value, with about twice the uncertainty.
Clearly the theoretical error on a HVP;LO μ in Eq. (5.1) is still too large to draw any conclusions regarding the presence of new physics and must be reduced by around a factor of 10 to reach the 0.2% target precision of the Muon g − 2 Experiment. Three key ingredients are still missing from our calculation of a HVP;LO μ described here: the effect of the difference between the u-and d-quark masses and of the quarks' electric charges on the light-quark connected contribution, and the contribution to the total from quark-disconnected diagrams involving u, d, s, and c quarks. Work on all of these is in progress [28,61]. Because they are all small corrections, however, relatively high accuracy is not needed. Ultimately calculations will be done on gluon-field configurations in which the sea quarks have both color and electric charges. Generation of such an ensemble is underway [62].
We must also further reduce the uncertainty on the lightquark connected contribution a ll μ ðconn:Þ in Eq. (5.4). The error budget (Table IV) is dominated by the lattice-spacing uncertainty, statistical errors, and the continuum extrapolation. The last two can be reduced by increasing statistics, so that the results at each lattice spacing value are more precise, and hence provide better constraints on the continuum extrapolation. We have demonstrated here that a calculation with nearly 0.5 million correlators (our highstatistics sample at a ¼ 0.15 fm) resolves issues around how to handle statistical uncertainties at large Euclidean times. Such a sample is numerically expensive to obtain on finer lattices, although tripling the statistics is certainly feasible using the truncated solver method. We estimate that this would reduce our total uncertainty to 1%. Further improvements may be achieved by analyzing additional correlation functions that include two-pion operators to better resolve the large-time behavior of the vector-current correlation functions [63,64]. To get below 1% requires a reduction in the uncertainty on the physical value of w 0 that determines the lattice spacing (w 0 =a is determined very precisely; see Table I). This uncertainty currently relies on a determination of the pion decay constant, f π , on the lattice [38]. The error budget in [38] shows that the dominant uncertainties are related to statistical precision and extrapolation to the physical point where w 0 f π is fixed against experiment (assuming a value of V ud from nuclear physics). An improvement by a factor of 2 in this uncertainty seems feasible with the higher statistics gluon-field ensembles now available with physical m u=d on finer lattices. Analysis on QCD þ QED gluon field ensembles will be important here too to take into account fully the fact that the decaying pion is an electrically charged particle. We also plan to investigate other quantities for determining the lattice spacing.
Given the above discussion, a reduction in uncertainty on the lattice-QCD result for the hadronic vacuum polarization contribution to the muon g − 2 to ≈0.5% is certainly feasible on the timescale of the new experiments. This would give precision comparable to that currently available from using experimental information on e þ e − → hadrons and would allow lattice-QCD results to play a significant role in the unfolding story of the search for new physics in the anomalous magnetic moment of the muon. FIG. 13. Comparison of our result in Eq. (5.1) for the leadingorder hadronic-vacuum-polarization contribution to the muon anomalous magnetic moment (magenta square) with results from N f ≥ 2 lattice QCD [14][15][16][17]23,48,54] (blue and purple squares) and from experimental e þ e − cross-section data [5][6][7]55] (red and orange triangles). The filled black circle shows the value of a HVP;LO μ that is implied by the measurement of a μ by BNL experiment E821 [1] assuming no contributions beyond the Standard Model; vertical dashed lines denote the AE1σ range [25].