Light quark vacuum polarization at the physical point and contribution to the muon g − 2

Santiago Peris Department of Physics and IFAE-BIST, Universitat Autònoma de Barcelona, Barcelona, E-08193 Bellaterra, Spain Abstract We report on the computation of the connected light quark vacuum polarization with 2+1+1 flavors of HISQ fermions at the physical point and its contribution to the muon anomalous magnetic moment. Three ensembles, generated by the MILC collaboration, are used to take the continuum limit. The finite volume correction to this result is computed in the (Euclidean) time-momentum representation to NNLO in chiral perturbation theory. We find all μ(HVP) = (659 ± 20 ± 5 ± 5 ± 4)× 10−10, where the errors are statistical and estimates of residual uncertainties from taking the continuum limit, scale setting, and truncation of chiral perturbation theory at NNLO. We compare our results with recent ones in the literature.


I. INTRODUCTION
Fermilab experiment E989 is measuring the anomalous magnetic moment of the muon (a µ = (g − 2)/2) with the goal of reducing the error on the BNL E821 [1] result by a factor of four. An upcoming experiment at J-PARC, E34, aims to do the same with a completely different technique. Lattice calculations of the hadronic contributions to the muon g −2, like the one reported here, are crucial to obtain and cross-check the Standard Model value to the same accuracy in order to discover new physics or lay to rest the longstanding discrepancy between theory and experiment.
In this paper we focus on the leading hadronic vacuum polarization (HVP) contribution to the muon anomaly. The aim is to test the efficacy of modern noise-reduction techniques to reduce the statistical errors of Monte Carlo methods used in lattice QCD in the context of the HVP and to provide accurate finite volume corrections to these results using chiral perturbation theory at two-loop order.
The total HVP contribution to a µ comes from both connected-and disconnected-quark line diagrams shown in Fig. 1, for each flavor of quark in Nature. The u, d quark-connected contributions are by far the largest, and we only compute them in this work. Comparison to other recent precise calculations [2][3][4][5][6][7] will provide important validation for the lattice method. The plan of this paper is the following. In Sec. II we review the theoretical framework for the calculation, including important details of the lattice calculation and the calculation in chiral perturbation theory in Euclidean space of the leading and next-to-leading finite volume corrections to the HVP contribution to the muon g − 2. Section III presents our results and comparison to other calculations. In Sec. IV we give a summary of this work and discuss implications for future work and the important upcoming comparison with experiment. The Appendix reports details of the NNLO chiral perturbation theory calculation.

II. THEORETICAL FRAMEWORK
Using lattice QCD and continuum, infinite-volume (perturbative) QED, one can calculate the hadronic vacuum polarization (HVP) contribution to the muon anomalous magnetic moment [8][9][10], m µ is the muon mass, andΠ(q 2 ) is the subtracted HVP,Π(q 2 ) = Π(q 2 ) − Π(0), computed directly on a Euclidean space-time lattice from the Fourier transform of the vector current two-point function, = Π(q 2 )(−q µ q ν + q 2 δ µν ), j µ (x) is the electromagnetic current, and Q i is the quark electric charge in units of the electron charge e (the sum is over active flavors). The form in the second equation is dictated by Lorentz and gauge symmetries.
In the following it is convenient to use the time-momentum representation [11] which results from interchanging the order of the Fourier transform and momentum integrals in Eqs. (4) and (1), respectively.
C(t) = 1 3 x,i where C(t) is the Euclidean time correlation function, averaged over spatial directions, and T is the temporal size of the lattice, and a HVP µ is obtained in the limit T → ∞. We have anticipated the use of the lattice with a discrete version of Eq. (10). The weight w(t) is sometimes modified by replacing the continuum Euclidean momentum-squared with its lattice version [3]:ŵ Note the double subtraction [11][12][13] in the cosine term in Eq. (7): t 2 /2 cancels Π(0) "configuration-by-configuration" while the leading finite size correction is killed by the "-1".
The latter arises since Π µν (q 2 ) does not vanish as q 2 → 0 when the time extent of the lattice is finite [11], but instead leads to a thermal electric susceptibility. In fact such terms are not constrained by the Ward-Takahashi Identity which in infinite volume leads to Eq. (5) and are allowed by the lattice symmetries [11,13].

A. Finite volume chiral perturbation theory
In this section, we consider the calculation of finite-volume effects in a HVP µ to two loops, or next-to-next-to-leading order (NNLO) in chiral perturbation theory (ChPT), with the aim of correcting our lattice result for a HVP µ for finite-volume effects. With our pion masses near the physical value, it is safe to assume that even at NNLO the most significant finitevolume correction will come from pion loops, and we can thus restrict our calculation to isospin-symmetric two-flavor ChPT.
There are two possible strategies for doing this. One is to first carry out a continuum extrapolation, and using results from continuum ChPT to correct for finite-volume effects.
The other is to correct the results at each lattice spacing, to obtain infinite-volume results at fixed lattice spacing. As we are using staggered fermions, the second strategy requires the use of staggered ChPT (SChPT) [14,15]. If all our ensembles were at the same pion mass and volume, the two methods should yield equivalent results. However, both the pion masses and volumes of the three ensembles are slightly different (cf. Table I). In this case, applying the finite-volume correction at a fixed lattice spacing has the advantage that this automatically corrects for the slightly different volumes. 1 While a full two-loop SChPT calculation is outside the scope of this paper, it is easy to change the NLO continuum ChPT result into a SChPT result; one only has to carry out a weighted average over the different taste pion masses for a given ensemble [13]. In practice, what we will do is to first correct the finite-volume lattice results for a HVP µ using NLO SChPT, then extrapolate to the continuum limit, after which we apply the remaining NNLO continuum ChPT correction.
Because of the slight mistunings of the pion masses and volumes, there will be a systematic error associated with this last step, but this systematic error will be much smaller than it would be if we were to extrapolate to the continuum first, and then apply NLO plus NNLO continuum ChPT to correct for finite-volume effects.
While the vacuum polarization in finite volume to two loops has been calculated before in momentum space [17], 2 we will directly calculate C(t), defined in Eq. (6), in the time-momentum representation, for t > 0, in a spatial volume of linear size L, with periodic boundary conditions. 3 This makes the calculation somewhat simpler, because we do not have to consider diagrams that lead to contributions proportional to δ(t) (which, in momentum space, correspond to contact terms). Our result will depend on only two lowenergy constants, F , the pion decay constant in the chiral limit, and 6 , which is an order-p 4 low-energy constant appearing in the EM current at this order. 4 Of course, the ChPT expression for C(t) is only reliable for large t, while C(t) for all t > 0 is needed in the sum (8). 5 However, as already observed in Ref. [13], finite-volume effects are a long-distance effect, and one thus expects the finite-volume correction to this correlation function to be reliably estimated for all t > 0, so that we can, in fact, estimate the finite-volume effect in a HVP µ using ChPT. An advantage is that this avoids using models to go beyond NLO ChPT (which is the same as scalar QED), as was proposed in Ref. [20].
As we will see, the ChPT result for the difference 1 But not the slightly different pion masses [16]. 2 For recent work on finite-volume effects of order exp[−m π L] not using ChPT, see Ref. [18]. 3 We take the time extent to be infinite. 4 We use the notation and conventions of Ref. [19] for low-energy constants. 5 C(0) is not needed as the weight w(t) ∝ t 4 for small t.
is indeed well defined. 6 in what follows.
The pion contribution to the EM current, to the order we need, is given by 7 Working in Euclidean space, a relatively straightforward calculation in the time-momentum representation yields the result for C(t) to NNLO as in which and the sums over p and k are over the momenta 2π n/L, n i integer, in a box with periodic boundary conditions. In Eq. (14) we gave the result in d = 3 + spatial dimensions in order to regulate the UV divergence present in d = 3. After defining a renormalized r 6 by the limit d → 3 can be taken, yielding a finite result for C(t). The factor 10/9 is needed to isolate the light-quark connected part [13,20,21]. The pion mass m π appearing in Eq. (14) is the renormalized (physical) pion mass. This renormalization absorbs the lowenergy constants 3,4 which appear in the explicit calculation. We note that the terms in the double sum on the second line of Eq. (14) with p 2 = k 2 lead to a term proportional to t e −2Ept , leading to the expected energy shift for two pions in an I = 1, = 1 state in a finite volume [22]. 6 In general we define ∆f (L) = lim L→∞ f (L) − f (L) 7 There are contributions from other order-p 4 low-energy constants, but they do not appear in the result for C(t) after mass renormalization.
In order to extract the finite-volume corrections, we use the Poisson resummation formula Let us work out the extraction of ∆a HVP µ to NLO in ChPT, relegating the treatment of the NNLO contribution to the appendix. The NLO part of C(t) is obtained by dropping all terms of order 1/F 2 in Eq. (14). Employing Eq. (18), the NLO part C NLO (t) can be written as C NLO (t) = − 10 9 where n 2 is summed over all non-negative integers and [22] Z 00 (0, The n 2 = 0 term, with sin (npL)/(nL) → p and Z 00 (0, 0) = −1, yields C NLO (t) in the infinite-volume limit. Inserting Eq. (19) into Eq. (10) with T = ∞ and replacing the sum over t by an integral, we find that ∆a HVP, NLO µ = 10 9 with Using the parameter values of Table I, we the NLO contribution. We then assume that the next order in ChPT, which we did not compute, is again of order 0.4-0.45 times the NNLO contribution, and we use this estimate as our error.
The fact that the three values in Eq. (24) are different is due to the mistuning of the pion masses and volumes of the three ensembles. If we were to apply the correction to the continuum extrapolated value of a HVP µ , we would thus have to use some average, and the spread of 5.1 × 10 −10 between the three values would represent a systematic error associated with the mistuning. If we were to apply only the NNLO correction in the continuum limit, that spread would be reduced to 1.7 × 10 −10 (cf. Eq. (52)). Hence, as explained above, what we will do is to first use NLO SChPT to correct the value of a HVP µ at each lattice spacing, then extrapolate, and finally apply the NNLO correction computed in Eq. (52) in the Appendix.
Using the taste-split pion spectrum for each ensemble, 8 Finally, the n 2 = 0 term in Eq. (19) gives us access to the effect of taste breaking in the pion masses in infinite volume, to NLO in ChPT. We use this to compute the corresponding corrections for each of our ensembles, finding these to be equal to These corrections are to be added to the lattice result to correct for taste breaking in the pion spectrum in infinite volume, to NLO in ChPT. Of course, since taste breaking is a latticespacing effect, whether one adds these corrections or not should not matter in the continuum limit. The difference one finds between values extrapolated to the continuum limit with or without this correction thus gives an estimate of the systematic error associated with taking the continuum limit. 8 We thank Doug Toussaint for providing the pion spectra, and for discussions of the taste splittings.
We take a moment to describe the low-mode structure of the staggered fermion Dirac operator which plays a central role. For valence quarks we use the HISQ fermion Dirac operator minus the Naik term, so the following, which is true in general for naive staggered fermions, applies here. The staggered Dirac operator is the sum of a hermitian mass term which commutes with an anti-hermitian hopping term, so it satisfies (using even-odd ordering where m is the quark mass and M oe hops quarks from even to odd sites. Similarly, the Eigenvectors of the preconditioned operator are eigenvectors of M with squared magnitude eigenvalue, and the even part can be obtained from the odd part, The eigenvalues come in ± pairs: If n + = (n o , n e ) is an eigenvector with eigenvalue λ n , then is also an eigenvector with eigenvalue −λ n : Thus we can construct pairs of eigenvectors, n + , n − , corresponding to ±iλ for each eigen-pair (λ 2 , n o ) computed with the Lanczos algorithm.
The full-volume LMA takes advantage of the spectral decomposition of the quark propagator that requires only two independent volume sums instead of a volume-squared sum in the correlation function. We employ a conserved current (again, minus the three-hop Naik term) which makes the "meson fields" a bit more complicated, χ(x) are single component staggered fermion fields whose spinor nature is encoded in the staggered phases, η(x), arising from the spin diagonalization of the fermion action. The gauge links U µ (x) ensure the point-split current is gauge invariant. A spectral decomposition of the low-mode part of the quark propagator is used in the AMA and LMA procedures, where N low is the number of low modes. The two point, current-current correlation function then becomes where λ n is shorthand for either m ± iλ n , and the sums over eigenvectors run up to 2N low .
To compute the above we construct meson fields, with eigenvector ordering λ 0 , −λ 0 , λ 1 , −λ 1 , . . . , λ N low , −λ N low . The factor (−1) (m+n)x+m arises from the construction of n − from n + since even m or n always corresponds to n + while odd corresponds to n − .
The AMA and LMA procedures are used to produce an improved estimator for the expectation value of any observable O by adding and subtracting terms that are exactly equal in the infinite statistics limit. Outside this limit the unimproved and improved estimates are statistically equivalent, with the latter having smaller errors. Here the gain in statistical precision is assumed to apply after using the same computational resource. The combined AMA and full-volume LMA improved estimator is given by The first three terms on the right hand side of Eq. (34) correspond to AMA [23,24] while the last two supplement this with the full-volume LMA [3,[25][26][27][28]. The expensive "exact" (to numerical precision) calculation is done relatively seldom while the inexpensive "approx" calculation is done often to reduce the statistical error. The difference of the first term with the 2nd and 4th terms corrects the bias induced by the 3rd and 5th approximate terms.
Note that in this work the first two sums in Eq. (34) [29]. "LM" is the number of low-modes of the precon- In Fig. 2 the summand in Eq. (10) for each ensemble is shown along with the full volume LMA and AMA contributions. In the figure "total" refers to the sum of five terms in Eq. (34). As observed in Ref. [3] there is a huge reduction in statistical error from the low-mode average, the last term in Eq. (34) (compare the total with full volume LMA and without, which is just AMA). The error reduction is especially large for large distance, as expected since the low-modes dominate this region. For the 96 3 ensemble the number of low modes used was 2000 (×2) compared to the other two ensembles (2 × 3000), due to computer and memory resource limitations. This is unfortunate as one can see from Fig. 2 that the full volume LMA is not as effective. Even though it appears that the low-mode contribution is mostly saturated (since it is comparable for all three ensembles), apparently the extra low-modes for the two coarser ensembles are very effective at reducing statistical noise.
In order to reduce further the statistical errors on the integrated result, we employ the bounding method [2,3] wherein C(t), for t > T , is given by C(t) = 0 (lower bound), and i.e, the lowest (two pion) energy state in the vector channel. At sufficiently large T the bounds overlap, and an estimate for a µ can be made which may be more precise than simply summing over the noisy long-distance tail. In Fig. 3   data must be corrected for finite volume effects, taste symmetry breaking, and pion mass mistunings (similar corrections were made in Ref. [5]). To make the various corrections we employ the following general procedure. The contribution to a HVP To explore a more precise comparison with other results, we adopt the window method of Ref. [3]: where t 1 − t 0 is the size of the window and ∆ is a suitably chosen width that smears out the window at either edge. We choose windows to avoid both lattice artifacts at short distance and large statistical errors at long distance. Results for several windows and both weighting functions are tabulated in Tab. IV.
In Fig. 6 several continuum limits are shown for the window with t 0 = 0.4, t 1 = 1, and ∆ = 0.15 fm. For this window the statistical errors for each ensemble are very small, so it allows a precise regime to explore and understand discretization effects. Here we also ignore mass re-tunings and finite volume effects because they have a negligible effect, with the   Squares (crosses) correspond to uncorrected data points with weighting function w (ŵ); filled circles are taste-breaking corrected to NLO of w data points. Solid curves show linear fits in a 2 ; all three agree very well in the continuum limit. Dashed curves denote a fully constrained parametrization (no degrees of freedom) using both a 2 and a 4 terms. be true). However, we do investigate taste-breaking effects since these are significant. The lower two curves in Fig. 6 correspond to uncorrected data points and weighting functions w andŵ. At non-zero lattice spacing there is a noticeable effect, but the continuum limits are the same (see the last row in Tab. IV). Including the taste-breaking corrections shifts the data further, essentially making the curve flat, but the continuum limit is barely affected.
We also show totally constrained "fits," including an a 4 term, which lower the continuum limit slightly while significantly increasing the statistical error. Figure 7 displays results for two representative windows along with values from the recent RBC/UKQCD computation using domain wall fermions 11 . The results should agree in the continuum limit up to small systematics. We also show the corresponding dispersive/e + e − value, using the R-ratio compilation of Ref. [32]. The HISQ results lie above the DWF and dispersive ones. The largest difference is about 7 × 10 −10 , or roughly one percent of the total HVP contribution to a µ . Adding an a 4 term or leaving out the largest lattice spacing point tend to give somewhat lower values with larger statistical errors. Given the uncertainties it is difficult to conclude there is a significant discrepancy, though the spread seems uncomfortably large. It is interesting to note that the HISQ and DWF lattice spacing errors are comparable. Finite volume errors, which have not been included, are very small in the 0.4-1.0 fm window. Likewise, the absence of charm sea quarks in the DWF result is estimated from perturbation theory to be very small [3]. A third, smaller, lattice spacing ensemble is being generated by the RBC/UKQCD collaborations [33], which could firmly establish whether or not a discrepancy exists. Other groups could apply the window method to their existing data which would also be helpful.
A final check included for completeness comes from moments of the correlation function [31], For the first moment we find 0.0797 (27), 0.0841(39), and 0.069(39) for the three different ensembles, coarsest to finest, respectively. A simple linear extrapolation in a 2 yields Π ll 1 = 0.0884(86) which is consistent with the values in Refs. [3,5].

IV. CONCLUSION
We have presented a lattice QCD computation of the light quark HVP contribution to the muon anomaly with 2+1+1 flavors of HISQ fermions. Three ensembles at the physical point, generated by the MILC Collaboration, were used to take the continuum limit at fixed volume (L ≈ 5.5 fm), and the results are broadly consistent with those in the literature.
Using the window method, a precise comparison yields values that are a bit higher than the dispersive result and a recent one using DWF. Given the statistical and systematic errors it is not clear that a real discrepancy exists: a decisive determination requires additional computations.
Overall the statistical errors in this study are at the larger end of the range from recent studies [2-5, 10, 34]. This is primarily due to the fewer number of low modes and measurements on the largest lattice used in our study. We note that the quoted statistical errors on the continuum limit value given in Ref. [5], which uses the same ensembles as used here, are roughly three times smaller than a naive linear fit in a 2 would yield due to the use of a prior constraint on the slope. We opted not to use a prior, which could similarly decrease our errors. Therefore we have ignored smaller systematic shifts due to isospin breaking [3,5,35].
Nevertheless the error reduction techniques used here are demonstrably powerful. Future computations with more measurements, and in particular, that use more low modes, can have an impact.
We have also presented a calculation in chiral perturbation theory, in Euclidean space, through NNLO of the finite volume corrections to the HVP contribution to the muon g − 2.
The NNLO correction is large (∼ 1%) for physical pion mass and the lattice sizes used in current calculations, so it must be included for a precise comparison to experiment.
The computations presented here are important for the test of the Standard Model against the ongoing experiment at Fermilab and an upcoming one at J-PARC.

V. ACKNOWLEDGMENTS
We thank Ruth Van de Water for discussions on the continuum limit fit in Ref. [5]. This The software used for this work includes CPS and Grid.
Appendix: NNLO finite-volume correction At NNLO, using the resummation formula (18), C(t) of Eq. (14) can be rewritten as We note that, despite the appearance of k 2 − p 2 in the denominator in various places, this is always accompanied by a numerator that vanishes at k 2 = p 2 , and all functions we integrate over k and p are continuous. An implication is that if we (as we will do below) break up some of the terms containing the factor 1/( k 2 − p 2 ), any contributions from the apparent pole at k 2 = p 2 should be dropped. We will always regulate such poles such that they do not contribute to the integrals.
The first two lines give the infinite-volume result, while the remaining lines represent finite-volume corrections. These finite-volume corrections can be rearranged as in which the renormalization-group invariant 6 is defined by and where the limit d → 3 has already been taken in the first term. The first term in Eq. The first term (first two lines) of Eq. (40) can be dealt with in the same way as the NLO contribution; all one needs to do is to insert the expression between square brackets inside the integral over p in Eq. (21). Numerically, using For the third and fourth lines in Eq. (40), we need the integral which converges for d < 3. We use Cauchy's theorem to rewrite the k integral as an integral along the discontinuity of the square root across the cut which we choose along the positive imaginary axis starting at +im π . 12 The result is finite in the limit d → 3, and Eq. (44) then becomes equal to Z 00 (0, n 2 ) nm π L K 1 (nm π L) .
The numerical value of this expression is equal to 0.00399, 0.00375 and 0.00282 for the 96 3 , Here we again closed the contour in the upper half k plane, and regulated the poles at k = ±p − iη such that they are located in the lower half k plane, and thus do not contribute; cf. the explanation below Eq. (39).