Higher-order hadronic-vacuum-polarization contribution to the muon g-2 from lattice QCD

We introduce a new method for calculating the ${\rm O}(\alpha^3)$ hadronic-vacuum-polarization contribution to the muon anomalous magnetic moment from ${ab-initio}$ lattice QCD. We first derive expressions suitable for computing the higher-order contributions either from the renormalized vacuum polarization function $\hat\Pi(q^2)$, or directly from the lattice vector-current correlator in Euclidean space. We then demonstrate the approach using previously-published results for the Taylor coefficients of $\hat\Pi(q^2)$ that were obtained on four-flavor QCD gauge-field configurations with physical light-quark masses. We obtain $10^{10} a_\mu^{\rm HVP,HO} = -9.3(1.3)$, in agreement with, but with a larger uncertainty than, determinations from $e^+e^- \to {\rm hadrons}$ data plus dispersion relations.

= −9.3(1.3), in agreement with determinations from e + e − → hadrons data plus dispersion relations. The total uncertainty is below the target precision of the Muon g − 2 Experiment. We also provide new expressions suitable for computing the O(α 3 ) hadronic vacuum polarization contributions from the renormalized vacuum polarization function Π(q 2 ), or directly from the lattice vector-current correlator in Euclidean space.

I. INTRODUCTION
The anomalous magnetic moment of the muon (g µ − 2) is one of the most precisely-determined observables in particle physics, having been measured with an uncertainty of 0.54 parts-per-million by BNL Experiment E821 [1]. Because of this high experimental precision, and because the anomaly is mediated by quantummechanical loops in the Standard Model, the muon g µ −2 provides stringent constraints on new heavy or weaklycoupled particles. The present Standard-Model theory value lies below the BNL E821 measurement by more than three standard deviations [2]. To identify definitively whether this deviation is due to new particles or forces, both the theory and measurement errors must be improved. The Muon g µ − 2 Experiment recently began running at Fermilab, and aims to reduce experimental error by a factor of four [3]. In parallel, numerous efforts are underway by the lattice-QCD community to tackle the Standard-Model hadronic contributions [4][5][6][7][8][9][10][11][12][13], which are the largest source of theory uncertainty [2].
The largest source of uncertainty in the Standard-Model g µ − 2 is from the O(α 2 ) hadronic vacuumpolarization (HVP) contribution [2], a HVP µ , which is shown in Fig. 1. 1 This contribution can be obtained by combining experimental measurements of electronpositron inclusive scattering into hadrons with dispersion relations, and recent determinations from this approach quote errors of 0.4-0.6% [14][15][16]. The most precise calculation of the leading-order a HVP µ to-date from Ref. [8] employed four-flavor lattice QCD with physical- * ruthv@fnal.gov 1 The symbol α always denotes the electromagnetic coupling in this work.
mass pions to achieve a total error of ∼ 2%. A significant source of systematic uncertainty in this and all lattice-QCD results to-date is from the use of degenerate up-and down-quark masses; phenomenological estimates of this error are about 1% [17][18][19]. Recently, we calculated the strong-isospin-breaking correction to the leading-order, light-quark-connected contribution to a HVP µ directly for the first time with the physical values of m u and m d , thereby removing this important uncertainty contribution [20]. To match the target experimental precision, however, the error on a HVP µ must be further reduced to about 0.2%.
The O(α 3 ) "higher-order" hadronic vacuumpolarization contribution to g µ − 2 is roughly 1.5% that of the leading-order HVP contribution [2], and therefore only needs to be determined to around 10% to match the projected experimental precision. Experimental determinations from combining electron-positron inclusive scattering into hadrons data with dispersion relations quote errors of 0.4-0.9% [14,16,21]. Nevertheless, it is important to check these phenomenological values with ab-inito QCD calculations. Moreover, if the disagreement between theory and experiment persists or grows with the new Muon g µ − 2 measurement, a complete first-principles Standard-Model theory value will be essential for drawing conclusions about the presence or nature of new physics.
In this paper we calculate the higher-order HVP contribution to a HVP µ for the first time in lattice QCD. To enable us to focus on the methodology and error analysis, we use previously published lattice-QCD results for the Taylor coefficients of the renormalized vacuum polarization function ( Π(Q 2 )) from Refs. [8,[22][23][24] to construct both Padé [23] and Mellin-Barnes approximants [25] for Π(Q 2 ). Details on the lattice-QCD calculations can be found in these works. This paper is organized as follows. In Sec. II, we provide theoretical background on the hadronic-vacuumpolarization contributions to g µ − 2, and discuss our method for calculating the higher-order contributions. Next, in Sec. III we present our analysis and error budget. Last, in Sec. IV, we show our final result for a HVP,HO µ and compare with non-lattice determinations. Appendix A provides expressions suitable for computing the O(α 3 ) hadronic vacuum-polarization contribution to a HVP µ directly from lattice-QCD simulations, while App. B provides the definition of the N = 2 + 1 + 1 Mellin-Barnes approximant for the Π(Q 2 ) used in this paper. For completeness, App. C gives the values of the quark-connected Taylor coefficients employed in our analysis.

II. THEORETICAL BACKGROUND
The leading hadronic contribution to the muon anomalous magnetic moment arises from QCD corrections to the internal photon propagator in the O(α 2 ) one-loop muon vertex diagram, as shown in Fig. 1. At O(α 3 ), higher-order hadronic contributions arise from adding a second internal photon line (as in Fig. 2 (a)), adding a lepton loop to the existing photon line (as in Figs. 2 (a) and (b)), or adding a second insertion of the hadronic vacuum polarization bubble on the photon line (as in Fig. 2 (c)). Both the leading-and NLO HVP contributions can be obtained, with the help of dispersion relations, from the energy scan of the experimental "R-ratio" [14][15][16]21]: where s is the square of the center-of-mass energy. Table I two recent evaluations of the leading contribution and the individual higher-order contributions from diagrams (a), (b), and (c) by Jegerlehner [14] and Keshavarzi et al. [16]. The higher-order contributions are roughly 1.5% of the leading contribution, and do not contribute substantially to the total error on the Standard-Model theory value for a µ . Integrals for the O(α 3 ) contributions from diagrams (a)-(c) have been presented in the literature in terms of R γ (s) [26,27]. These formulations, however, are not suited for our use, particularly in the case of contribution (a). We therefore provide in Appendix A new expressions for these contributions that are amenable to use with lattice-QCD data. For each contribution, we provide two formulations to obtain a (i) First, we use the following relationship between R γ (s) and the renormalized vacuum polarization function [28], to derive expressions in terms of the renormalized vacuum polarization function Π(Q 2 ) ≡ Π(Q 2 ) − Π(0). 2 These are the higher-order analogs of the original Blum formula for the leading HVP contribution [29], and are given in Eqs. (A1), (A6), and (A13). We also provide expressions for the contributions from diagrams (a)-(c) directly in terms of the Euclidean vector-current correlator at zero momentum G(t) using the relationship between Π(Q 2 ) and G(t) below [28]: These are the higher-order analogs of the timemomentum representation formulated by Bernecker and Meyer for the leading HVP contribution, and are given in Eqs. (A3), (A11), and (A14). The higher-order HVP contributions are sensitive to the value of the renormalized vacuum polarization function at larger values of Q 2 than the leading-order contribution. Figure 3, left, plots the integrands for the leading-order and higher-order contributions as a function of Q 2 using the N = 2 + 1 + 1 Mellin-Barnes approximant for Π(Q 2 ) from Ref. [25]. The integrand for the leading-order contribution is also shown for comparison. The integrand of contribution (a) has large positive and negative contributions below Q 2 = m 2 µ that cancel substantially. Because of this, the large-Q 2 region is numerically important, with about 5% of the value of a it decreases less rapidly with Q 2 than the other contributions; about 10% of the value of a (c) µ comes from Q 2 > 1GeV 2 . Thus, it is important to employ approximants of Π(Q 2 ) that accurately reproduce the large-Q 2 behavior when calculating the higher-order contributions to a HVP µ .
The higher-order HVP contributions are sensitive to the value of the Euclidean-time correlator at similar times as the leading-order contribution. Figure 3, right, plots the integrands for the leading-order and higher-order contributions (a) and (b) as a function of correlator time t using G(t) obtained from the spectral representation of R γ (s). (The kernel for contribution (c) depends upon the product of the correlator at two times G(t)G(t ) and thus the integrand cannot be conveyed in a one-dimensional plot.) The leading-order (higher-order) kernels are proportional to t (t 2 ) at small Euclidean times, and are proportional to 1/t (approach a constant )at large times, and the integrands all peak at around t ∼ 0.8-1.0 fm. The contributions to a HVP µ from correlator data beyond 4 fm, which is approximately half the temporal extent (or less) of lattices employed in recent g − 2 calculations, are about 0.5% or less [8,11,12,30].

III. ANALYSIS
In this section we calculate the O(α 3 ) contributions to a HVP µ from the diagrams in Fig 2. First, in Sec. III A, we describe the approximants of the renormalized vacuum function used to calculate the higher-order HVP contributions. Next, we calculate the quark-connected contribution from light and heavy quarks in Sec. III B Last, in Sec. III C, we estimate the size of the quark-disconnected contribution.

A. Approximants of Π(Q 2 )
We calculate the higher-order contributions to a HVP µ using both Padé and Mellin-Barnes approximants of the renormalized vacuum polarization function in the QED integrals given in Appendix A. Both approaches employ the Taylor coefficients Π i of Π(Q 2 ) expanded about Q 2 = 0: As observed in Ref. [23], the Π i are proportional to the time-moments of the vector-current correlation function, and can be computed with small statistical errors in lattice QCD. Further, with both the Padé and Mellin-Barnes approches, only the first few Taylor coefficients are needed to obtain the leading-order HVP with a subpercent systematic uncertainty associated with the parameterization of Π(Q 2 ) [8,25]. Following the method introduced by the HPQCD Collaboration [23], we construct the [n, m] Padé approximants for the renormalized hadronic vacuum polarization function from the Π i 's. The true result for Π(Q 2 ) is guaranteed to lie between the [n, n] and [n, n − 1] Padé approximants. For the leading-order HVP contribution, the Padé approximants provide a sufficiently accurate approximation of Π(Q 2 ) both at low and high  [25], which employs preliminary moments of Rγ(s) provided by Keshavarzi et al. [16]. The leading-order integrand is also shown as a solid magenta line for comparison. Right: integrands of Eqs. (A3) (blue squares) and (A11) (green diamonds) obtained from the parameterization of Rγ(s) provided by Jegerlehner in his public alphaQED FORTRAN package [31]. The leading-order integrand is also shown as magenta circles for comparison.
Q 2 that the associated uncertainty in a HVP µ is below 1% by n = 2 [8]. Unfortunately, however, one cannot use the [n, n − 1] approximants Π(Q 2 ) to calculate the contributions to a HVP µ from diagrams (a) and (c). This is because Π [n,n−1] (Q 2 ) ∼ Q 2 as Q 2 → ∞, making the integrals diverge in this limit. The integrals using the [n, n] Padé approximants are well behaved, but another approach is needed to quantify the uncertainty in the higher-order contributions to a HVP µ from the parameterization of Π(Q 2 ).
Recently de Rafael and Charles et al. introduced the method of "Mellin-Barnes approximants" to obtain a HVP µ from the Taylor coefficients of Π(Q 2 ) [25,32]. This approach uses the fact that the hadronic spectral function Im Π(q 2 )/π in QCD is positive and approaches a constant as Q 2 → ∞ to identify a class of functions that can be employed as successive approximants to the Mellin transform M(s) of the hadronic spectral function. Given N moments of the Mellin transform M(−n), the Mellin-Barnes approximant M N smoothly interpolates between these known values, and approaches the asymptotic value of M(s) from leading-order perturbative QCD as s → ∞. The Mellin moments are trivially related to the Taylor coefficients of Π(Q 2 ) as 2) The first term in the moment expansion of the hadronic spectral function provides a rigorous upper bound on Π(Q 2 ) and a HVP µ [33]. In practice, the N = 1 approximant obtained using M(0) from experimental R γ data yields a value for the leading-order HVP contribution that already agrees with the full result to better than 1% [25]. Figure 4 plots the Padé and Mellin-Barnes approximants for Π(Q 2 ) calculated from the first four moments of R γ (s) [16], and compares them with the exact result obtained from direct integration of R γ (s). The Mellin-Barnes approximants are closer to the exact Π(Q 2 ) than the Padés because they are constrained to satisfy the asymptotic perturbative-QCD behavior as Q 2 → ∞. Note, however, that the rate at which the Mellin-Barnes approximants approach the true Π(Q 2 ) depends upon the specific functional form employed at each order. In particular, the difference between successive approximants is not guaranteed to decrease with increasing N .

B. Quark-connected contribution
We calculate the O(α 3 ) quark-connected contribution to a HVP µ using the Taylor coefficients of Π(Q 2 ) obtained by the HPQCD Collaboration in Refs. [8,[22][23][24]. The u, d, and s-quark Taylor coefficients were calculated on the MILC Collaboration's QCD four-flavor gauge-field configurations with highly-improved staggered (HISQ) sea and valence quarks [35,36]. The b-quark Taylor coefficients were also calculated on the HISQ ensembles, but with a radiatively-improved nonrelativistic QCD action for the b quarks [37,38]. The c-quark Taylor coefficients were calculated with HISQ valence quarks, but on MILC's three-flavor ensembles with asqtad sea quarks [39][40][41]. The MILC ensembles are isospinsymmetric, i.e. the up and down sea-quark masses are  [16] and quoted in Ref. [25]. The exact result is shown as a solid black line for comparison [34].
degenerate. The light-quark mass varies from m l = m s /5 to Nature's value m l ∼ m s /27, making a chiral extrapolation unnecessary, and the strange-(and charm-) seaquark masses are fixed to close to their physical values. We employ light-and strange-quark Taylor coefficients on two ensembles with physical light-quark masses and lattice spacings a ≈ 0.15 fm and 0.12 fm from Refs. [8,23]. Table IV gives the light-and strange-quark connected Taylor coefficients used in our analysis. The values of Π (ud) i include corrections for the finite lattice spatial volume and and nonzero lattice spacing computed at one-pion-loop order within scalar QED [18]. We employ charm-and bottom-quark Taylor coefficients from Refs. [22,24], which provide values of Π  Table V gives the heavy-quark connected Taylor coefficients used in our analysis.
To calculate the connected contribution to a (HO) µ , we first sum the individual Taylor coefficients Π i , and then use the total to construct the Padé and Mellin-Barnes approximants for Π(Q 2 ). Beyond N = 2, the functional forms of the Mellin-Barnes approximants are not unique; Appendix B gives the form of Π 2+1+1 (Q 2 ) used here. We then use the resulting approximants for Π(Q 2 ) in the QED integrals, Eqs. (A1), (A6), and (A13), to obtain the quark-connected contributions to a HVP µ from the diagrams in Fig. 2. On each ensemble, and for each contribution (a)-(c), we average the values from the Padé and Mellin-Barnes approximants, and take half the difference between the two as the systematic uncertainty from the parameterization of Π(Q 2 ). ) hadronic-vacuum-polarization contributions to gµ −2 on two physical-mass HISQ ensembles obtained using [2,2] Padé and N = 2 + 1 + 1 Mellin-Barnes approximants for Π(Q 2 ). The uncertainties are from the errors on the Taylor coefficients and, for the averages, from the use of approximants for Π(Q 2 ).  Table II gives the results on the two ensembles employed in our analysis. Figure 5 shows the total O(α 3 ) quark-connected contribution to a HVP µ -obtained by summing contributions (a)-(c) in the rows labeled "average" in Table II -versus squared lattice spacing. The data do not display any significant lattice-spacing dependence, so we fit them to constant to obtain the continuum-limit value of a HVP,HO µ . We also consider an alternative linear extrapolation in a 2 to a function of the form a HVP,HO with Λ = 500 GeV a typical QCD scale. The linear-fit result for c a 2 is consistent with zero, and for a HVP,HO µ is close to the value from the constant fit. We therefore conclude that discretization effects are smaller than the fit error on a HVP,HO µ , and do not assign a separate systematic error from this source.
The HPQCD Collaboration reduced the statistical errors in the light-quark connected Taylor coefficients in Ref. [8] by using fit results for the vector-current correlators for times greater than 1.5 fm. Although the lowestenergy states in these correlators are I = 1 ππ pairs, no evidence of such states was seen in the two-point fits, and the ground-state energies obtained are consistent with the experimental ρ 0 meson mass. HPQCD estimate the contribution to the leading-order light-quark connected contribution to a HVP µ from the omitted ππ states within scalar QED to be 3 × 10 −10 . We expect ππ contributions to be similar in size for the dominant higher-order diagrams (a) and (b) because the integrands in Eqs. (A1) and (A6) are proportional to Π(Q 2 ), just as for the leading-order hadronic vacuum polarization. Hence, we take the same percentage error of 0.5% as the uncertainty in a HVP,HO laborations calculated the strong-isospin-breaking correction to a HVP µ for the first time with physical values m u and m d [20]. They obtain +1.5(7)% for the relative correction that should be applied to the leading-order lightquark connected contribution, in agreement with phenomenological estimates [17][18][19]. Here we use +1.5(1.0)% to correct the continuum-limit value of a HVP,HO µ from Fig. 5, where we have taken a larger uncertainty of 1% on the relative correction to account for the fact that the shift was not calculated directly for the higher-order hadronic vacuum polarization.
The QCD gauge-field ensembles employed in our analysis do not include effects due to the quarks' nonzero electromagnetic charges in Nature. The dominant QED effect in a HVP µ arises from producing a hadron polarization bubble consisting of a π 0 -γ pair. Following Hagiwara et al. [42] we calculate the contribution to a HVP,HO µ from e + e − → π 0 γ in the region 0.6 < √ s < 1.03 GeV using the latest experimental data for this channel from the SND Experiment [43]. We obtain which is approximately 0.6% of the total quark connected contribution. We therefore take 1% as the error from the omission of electromagnetism in the simulations. Finally, as discussed in Appendix A, in order to express higher-order contribution 2(a) in Fig. 1 in terms of the renormalized vacuum polarization function, we must drop terms in the original integrand [26,27] that are proportional to (m 2 µ /s) n log 2 (m 2 µ /s). We have calculated the numerical size of these terms from experimental R γ (s) data [31] and, although they are small, they are not negligible given the size of our statistical and other systematic uncertainties. To account for the omission of the "log 2 " in our calculation of contribution 2(a) via Eq. (A1), we therefore include an additional systematic uncertainty of 1 × 10 −10 , which is almost twice the size of these terms calculated from R γ (s) data. Table III gives the complete error budget for the O(α 3 ) quark-connected contribution to a HVP µ . The largest uncertainties are associated with the omitted "log 2 " terms in contribution 2(a) and from the use of Padé and Mellin-Barnes approximants for the renormalized vacuum polarization function. Although the estimated uncertainties from the omission of QED and isospin breaking in the gauge-field configurations, and from low-lying ππ states in the vector-current correlators, are based on calculations for the leading-order vacuum polarization, they are about four times smaller, and do not contribute substantially to the total error. We obtain for the quark- (3.5) where "lat." denotes the sum of contributions associated with the underlying lattice-QCD calculations of the Taylor coefficients.

C. Quark-disconnected contribution
Although several lattice-QCD calculations of the leading-order quark-disconnected contribution to a HVP µ are available [7,12,44], these publications do not provide the Taylor coefficients of the renormalized vacuum polarization function. 3 We therefore estimate the values of the quark-disconnected Taylor coefficients assuming ground-state dominance of the vector-current correlators 3 In Ref. [10], the BMW Collaboration provides the first two Taylor coefficients Π , which are not sufficient to construct the [2,2] Padé and N = 2 + 1 + 1 Mellin-Barnes approximants.
as in Ref. [44]. Using Eq. (11) of that work, (25), 0.78265(12)} GeV from the PDG [45] and {f ρ , f ω } = {0.21(1), 0.20(1)} GeV yields and similar results for the higher Taylor coefficients. Both the leading O(α 2 ) contribution to a HVP µ and the domiant O(α 3 ) contributions from diagrams (a) and (b) are proportional to the Taylor coefficient Π 1 at lowest order in the small-Q 2 expansion. Further, the dominant quark-connected contribution is from the light up and down quarks. We therefore take −1.3(1.2)% as the correction and uncertainty due to the omission of quark-disconnected contributions in our analysis. We note that our estimate in Eq (3.7) is consistent with recent lattice-QCD calculations of the leading-order quarkdisconnected contribution with physical-mass pions from the BMW [12] and RBC/UKQCD Collaborations [7], who obtain for the ratio a (LO,disc.) µ /a (LO,u/d conn.) µ approximately -2.0% and -1.5%, respectively.

IV. RESULT AND OUTLOOK
We obtain the total O(α 3 EM ) hadronic vacuum polarization contribution to g µ − 2 by adding our calculation of the the quark-connected contribution, Eq. where the first two errors errors are from the quarkconnected and quark-disconnected contributions, respectively. We list the error from omission of the "log 2 " terms separately, since it does not arise from the use of lattice QCD to obtain the renormalized vacuum polarization function. This error could be eliminated with a different trick for expressing contribution (a) in terms ofΠ(Q 2 ) than the one employed here. The uncertainty on the quark-connected contribution stems primarily from our use of Padé and Mellin-Barnes approximants forΠ(Q 2 ), which we employ so that we can exploit already-published values of the Taylor coefficients. We anticipate reducing this error in a future paper that also includes an update of our determination of the leading-order hadronic vacuum polarization contribution [8] by calculating the higherorder contributions directly from the lattice correlation functions using the alternative formulae in Appendix A, and by analyzing a larger data set with more ensembles and finer lattice spacings.
Our result in Eq. (4.1) is the first lattice-QCD determination of the higher-order hadronic vacuum polarization contribution to g µ − 2, and is consistent with determinations from e + e − → hadrons data [14,16,21]. Although the lattice-QCD uncertainty is approximately ten times larger than from experiment plus dispersion relations, it is still below the target uncertainty of the Muon g − 2 Experiment of 0.14 ppm, or δa µ ∼ 1.6 × 10 −10 [3]. It therefore provides a new necessary ingredient in reaching the goal of obtaining a purely ab-initio-QCD determination of the hadronic contributions to g µ − 2. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

Contribution (a)
A complete analytical result for the contribution from the diagrams in (a) of Fig. 2 was first presented by Barbieri and Remiddi in Ref. [26]; in this work they also provide an expansion to first order in m 2 µ /s. Later, in Ref. [27], Krause derived an asymptotic expansion for the kernel function in terms of the parameter r = m 2 µ /s, which is more amenable to numerical integration. We start with the asymptotic expression given in Eq. (7) of Krause, which contains powers and logarithms of r. Equation (7) does not have the form needed to exploit the relationship between R γ (s) and the renormalized vacuum polarization function in Eq. (2.2). As suggested by Groote et al. [46], however, one can exploit generating integral representations of r n and r n log(r) to express the pure polynomial and log terms in the asymptotic expansion of the kernel function in terms of Π. Using Eqs. (39)- (42) of that work, and discarding terms proportional to log 2 (r) yields the following integral expression for contribution (a) in terms of the renormalized vacuum polariza-tion function: with Checking the size of the omitted logarithmic terms using experimental data for R γ (s) [31], we find that they are below 1 × 10 −10 .
Alternatively, contribution (a) is given in terms of the Euclidean zero-momentum correlator by with and The factors of t 2 and 1/t 2 in Eqs. (A3) and (A4), respectively, are chosen to make the kernel functionK (a) (t) dimensionless. With these formulae, contribution (a) can be obtained from a simple weighted sum of G(t) as in the leading-order case.

Contribution (b)
We start from Eq. (9) of Ref. [27] and make the change of variables Q 2 = m 2 µ x 2 /(1 − x). The contribution from diagram (b) in Fig. 2 is then given in terms of the renormalized vacuum polarization function by where the lepton loop function is and K E (Q 2 ) is the standard kernel function introduced by Blum in Ref. [29]: Thus, the expression in Eq. (A6) is simply the leadingorder QED integral with the replacementΠ(Q 2 ) → 8πα × Π Q 2 F m 2 e , Q 2 . The analogous contribution from the τ lepton is negligible because it is suppressed by m 2 µ /m 2 τ .
Contribution (b) can also be obtained from a weighted sum of the Euclidean zero-momentum correlator as in the leading-order case: with the dimensionless kernel

Contribution (c)
We start from Eq. (13) of Ref. [27]. Diagram (c) in Fig. 2 contains two hadronic insertions, and thus the contribution depends upon the square of the renormalized vacuum polarization function: In this case, the expression in Eq. (A13) has the form of the 1-loop QED integral, but with the replacement Π(Q 2 ) → 4πα × Π(Q 2 ) 2 . When contribution (c) is expressed in terms of the Euclidean zero-momentum correlator, the two powers of the vacuum polarization function above yield two integrals over times t and t : with the dimensionless kernel This formulation is slower to implement numerically than the analogous formulae for contributions (a) and (b) due to the double integral.  [8,23]. The quoted errors include statistics, the uncertainty on the vector-current renormalization factor, the (correlated) uncertainty from setting the lattice spacing, and the uncertainty on the corrections. The factor of the quarks' electromagnetic charges (Q 2 u + Q 2 d ) is included in the definition of the Πjs.