The hadronic light-by-light scattering contribution to the muon anomalous magnetic moment from lattice QCD

2RIKEN BNL Research Center, Brookhaven National Laboratory, Upton, New York 11973, USA 3Physics Department, Columbia University, New York, New York 10027, USA 4Department of Physics, Nagoya University, Nagoya 464-8602, Japan 5Nishina Center, RIKEN, Wako, Saitama 351-0198, Japan 6Physics Department, Brookhaven National Laboratory, Upton, New York 11973, USA 7Universität Regensburg, Fakultät für Physik, 93040, Regensburg, Germany


I. INTRODUCTION
The anomalous magnetic moment of the muon is providing an important test of the Standard Model. The current discrepancy between experiment and theory stands between three and four standard deviations. An ongoing experiment at Fermilab (E989) and one planned at J-PARC (E34) aim to reduce the uncertainty of the BNL E821 value [1] by a factor of four, and similar efforts are underway on the theory side . A key part of the latter is to compute the hadronic light-by-light (HLbL) contribution from first principles using lattice QCD [31][32][33][34][35][36][37]. Such a calculation, with all errors under control, is crucial to interpret the anticipated improved experimental results.
The magnetic moment is an intrinsic property of a spin-1/2 particle, and is defined through its interaction with an external magnetic field B, H int = −µ · B. Here where S is the particle's spin, q and m are the electric charge and mass, respectively, and g is the Landé g-factor. The Dirac equation predicts that g = 2, exactly, so any difference from 2 must arise from interactions. Lorentz and gauge symmetries tightly constrain * ljin.luchang@gmail.com the form of the interactions, where J ν is the electromagnetic current, and F 1 and F 2 are form factors, giving the charge and magnetic moment at zero momentum transfer (q = p − p = 0), or static limit. u(p) andū(p) are Dirac spinors. The anomalous part of the magnetic moment is given by F 2 (0) alone, and is known as the anomaly, a µ ≡ (g − 2)/2 = F 2 (0).
The desired matrix element in (2) is extracted in quantum field theory from a correlation function of fields as depicted in the Feynman diagrams shown in Fig. 1. Here we work in coordinate (Euclidean) space and use Lattice QCD for the hadronic part which is intrinsically nonperturbative. QED is treated using the same discrete, finite, lattice as used for the hadronic part, while we remove the spatial zero modes of the photon propagator. This method is called QED L [38]. It is perturbative with respect to QED, i.e, only diagrams where the hadronic part is connected to the muon by three photons enter the calculation.

II. QEDL METHOD
Here the muon, photons, quarks, and gluons are treated on a finite, discrete lattice. The method is described in arXiv:1911.08123v1 [hep-lat] 19 Nov 2019 FIG. 1. Leading contributions from hadronic light-by-light scattering to the muon anomaly. The shaded circles represent quark loops containing QCD interactions to all orders. Horizontal lines represent muons. Quark-connected (left) and disconnected (right) diagrams are shown. Ellipsis denote diagrams obtained by permuting the photon contractions with the muons and diagrams with three and four quark loops with photon couplings (See Fig. 3).
detail in Ref. [32], and the diagrams to be computed are shown in Figs. 2 and 3. It is still not possible to do all of the sums over coordinate space vertices exactly with currently available compute resources. Therefore we resort to a hybrid method where two of the vertices on the hadronic loop(s) are summed stochastically: point source propagators from coordinates x and y are computed, and their sink points are contracted at the third internal vertex z and the external vertex x op . Since the propagators are calculated to all sink points, z and x op can be summed over the entire volume. The sums over vertices x and y are then done stochastically by computing many (O(1000)) random pairs of point source propagators. To do the sampling efficiently, the pairs are chosen with an empirical distribution designed to frequently probe the summand where it is large, less frequently where it is small. Since QCD has a mass-gap, we know the hadronic loop is exponentially suppressed according to the distance between any pair of vertices, including |x − y|. As we will see, the main contribution comes from distances less than about 1 fm. The muon line and photons are computed efficiently using FFT's; however, because they must be calculated many times, the cost is not negligible.
Two additional, but related, parts of the method bear mentioning. First, the form dictated by the right hand side of Eq. 2 suggests the limit q → 0 is unhelpful since the desired F 2 term is multiplied by 0. Second, in our Monte Carlo lattice QCD calculation the error on the F 2 contribution blows up in this limit. The former is avoided by evaluating the first moment with respect to x op at the external vertex and noticing that an induced extra term vanishes exponentially in the infinite volume limit [32]. This moment method allows the direct calculation of the correlation function at q = 0, and hence F 2 (0). To deal with the second issue, we first recall that it is the Ward identity that guarantees the unwanted term to vanish in the moment method. We thus enforce the Ward identity exactly on a configuration-by-configuration basis [32]. i.e., before averaging over gauge fields by inserting the external photon at all possible locations on the quark loop in Fig. 2. This makes the factor of q in Eq. (2) exact for each measurement and not just in the average and reduces the error on F 2 (0) significantly. Implementing the above techniques produces an order O(1000) fold • Point source photons at x and y.
• Importance sampling is used in choosing x and y.
− Major contribution comes from the region where x and y are not far separated.
− In fact, we can evaluate all possible (upto discrete symmetries) relative positions for distance less than a certain value rmax, which is normally set to be 5 lattice units.
• Moment method for xop. Evaluate F2(q 2 ) at q = 0 directly.  improvement in the statistical error over the original nonperturbative QED method used to compute the hadronic light-by-light scattering contribution [31].
The quark-disconnected diagrams that occur at O(α 3 ) are shown Fig. 3. All but the upper-leftmost diagram vanish in the SU (3) flavor limit and are suppressed by powers of m u,d − m s , depending on the number of quark loops with a single photon attached. For now we ignore them and concentrate on the leading disconnected diagram which is computed with a method [33] similar to the one described in the previous section. To ensure the loops are connected by gluons, explicit vacuum subtraction is required. However, in the leading diagram the moment at x op implies the left-hand loop in Fig. 3 vanishes due to parity symmetry, and the vacuum subtraction is done to reduce noise.
As for the connected case, two point sources (at y and z in Fig. 3) are chosen randomly, and the sink points are summed over. M propagators are computed, and all M 2 combinations are used to perform the stochastic sum. This "M 2 trick" is crucial to bring the statistical fluctuations of the disconnected diagram under control.

III. LATTICE SETUP
The simulation parameters are given in Tab. I. All particles have their physical masses (isospin breaking for the up and down quark masses is not included). The discrete Dirac operator is known as the (Möbius) domain wall fermion ((M)DWF)) operator. Similarly the discrete gluon action is given by the plaquette plus rectangle Iwasaki gauge action. Additionally, three ensembles with larger lattice spacing employ the dislocation-suppressingdeterminant-ratio (DSDR) to soften explicit chiral symmetry breaking effects for MDWFs [39]. We use All Mode Averaging (AMA) [44] and Multi-grid Lanczos [45] techniques to speed up the fermion propagator generation.
The muons and photons take discrete free-field forms.
The muons are DWFs with infinite size in the extra fifth dimension, and the photons are non-compact in the Feynman gauge. In the latter all modes with q = 0 are dropped, a finite volume formulation of QED known as QED L [38].  [40]. The lattice spacing a, spatial extent L, extra fifth dimension size Ls, muon pion mass mπ, and number of QCD configuration used for the connected and the disconnected diagrams.

IV. RESULTS
Before moving to the hadronic case, the method was tested in pure QED [32]. Results for several lattice spacings and box sizes are shown in Fig. 4. The systematic uncertainties are large, but under control. Note that the finite volume errors are polynomial in 1/L and not exponential, due to the photons which interact over a long range. The data are well fit to the form 4 .
The continuum and infinite volume limit is F 2 (0) = 46.9(2) stat × 10 −10 for the case where the lepton mass in the loop is the same as the muon mass, which is quite  48I con 64I con 24D con 32D con 48D con 32Dfine con 48I discon 64I discon 24D discon 32D discon 32Dfine discon FIG. 5. Cumulative contributions to the muon anomaly, connected (upper) and disconnected (lower). r is the distance between the two sampled currents in the hadronic loop (the other two currents are summed exactly) and the horizontal axis is the cumulative contributions from r and below. 24 3 IDSDR (squares), 24 3 IDSDR (squares), 32 3 IDSDR (crosses), 48 3 Iwasaki (diamonds), and 64 3 Iwasaki (plusses).
Our physical point calculation [33] started on the 48 3 , a −1 = 1.730 GeV, Iwasaki ensemble listed in the first column of Tab. I, for which we found a con µ = 11.60(0.96) stat × 10 −10 , a discon µ = −6.25(0.80) stat × 10 −10 , and a tot µ = 5.35(1.35) stat × 10 −10 for the connected, leading disconnected, and total HLbL contributions to the muon anomaly, respectively. The errors quoted are purely statistical. We have since improved the statistics on the leading disconnected diagram with measurements on 39 additional configurations, and the contribution becomes −6.03(60) × 10 −10 . Since then we have computed on several additional ensembles in order to take the continuum and infinite volume limits (see Tab. I).
The results are displayed in Fig. 6 along with curves obtained with the following equation: where a I , a D represent the lattice spacings for the Iwasaki and I-DSDR ensembles respectively. For the Iwasaki ensembles, we define the variable a D to be zero and vice versa. Therefore the lattice spacing is always equal to a = a I + a D . We allow different a 2 coefficients for the Iwasaki and I-DSDR ensembles as the gauge actions are different. The lattice spacings for the I-DSDR ensembles are not small enough to allow us to ignore the a 4 effects, and therefore we include them in the fit. As we only have two lattice spacings for the I-DSDR ensembles, with both a 2 and a 4 effects unknown, we cannot extrapolate to the continuum just with the I-DSDR ensembles. Therefore, based on this fit form, the continuum limit is obtained from the two Iwasaki ensembles, and the I-DSDR ensembles are used to obtain the volume dependence only. In particular, the 32Dfine ensemble does not affect the fitted a µ at all. It only helps to determine the parameter c D 2 , which provides evidence for the size of the potential O(a 4 ) systematic errors. We find for the connected, disconnected, and total contributions, a con µ = 23.76(3.96) stat (4.89) sys × 10 −10 , a discon µ = −17.12(3.46) stat (4.41) sys × 10 −10 , a tot µ = 6.80(4.65) stat (1.56) sys × 10 −10 , respectively. For the total contribution, we fit the total contribution for each ensemble, which is slightly different from the sum of the fitted results from the connected and the disconnected parts. Notice there is a large cancellation between the connected and disconnected diagrams that persists for a → 0 and L → ∞, so even though the individual contributions are relatively well resolved, the total is not. The cancellation is expected since hadronic light-by-light scattering at long distance is dominated by the π 0 which contributes to both diagrams, but with opposite sign [35,42,43]. Notice also that the a 2 and 1/L 2 corrections are individually large but also tend to cancel in the sum.
The systematic errors mostly result from the higher order discretization and finite volume effects which are not included in the fitting formula Eq. (5). We therefore estimate the errors through the change of the results after adding a corresponding term in the fitting formula. For O(1/L 3 ), we add another 1/(m µ L) 3 term with the same coefficient as the 1/(m µ L) 2 term. For O(a 4 ) effects, we add an a 4 term also for the Iwasaki ensembles with coefficient similar to the I-DSDR ensembles. For O(a 2 log(a 2 )) effects, we multiply the discretization effect terms in Eq. (5) by (1 − (α S /π) log(a 2 GeV)). For O(a 2 /L), we multiply the discretization effect terms in Eq. (5) by (1 − 1/(m µ L)). In addition, for the only two contributions which we have not included in the present HLbL calculation: (a) strange quark contribution to the connected diagrams; (b) sub-leading disconnected diagrams' contribution. We have performed lattice calculations with the QED ∞ approach [47] on the 24D ensemble to estimate the systematic errors. These systematic errors are added in quadrature and summarized in Tab. II. In the supplementary materials, these systematic errors are discussed in more detail.
While the large relative error on the total is a bit unsatisfactory, we emphasize that our result represents an important estimate on the hadronic light-by-light scattering contribution to the muon anomaly, with all systematic errors controlled. It appears that this contribution cannot bring the Standard Model and the E821 experiment in agreement.   6. Infinite volume extrapolation. Connected (top), disconnected (middle), and total (bottom). We have use the hybrid method to calculate the continuum limit for the connected contribution.
In fact we can do even a bit better with the data on hand. As seen in Fig. 5, which shows the cumulative sum of all contributions up to a given separation of the two sampled currents in the hadronic loop, the total connected contribution saturates at a distance of about 1 fm for all ensembles. This suggests the region r > ∼ 1 fm adds mostly noise and little signal, and the situation gets worse in the limits. A more accurate estimate can be obtained by taking the continuum limit for the sum up to r = 1 fm, and above that by taking the contribution from the relatively precise 48 3 ensemble. We include a systematic error on this long distance part since it is not extrapolated to a = 0. The infinite volume limit is taken as before. This hybrid procedure yields a con µ = 24.16(2.30) stat (5.12) sys × 10 −10 , with a statistical error that is roughly 2× smaller and the additional O(a 2 ) systematic error from the hybrid procedure is only 0.20 × 10 −10 . Unfortunately a similar procedure for the disconnected diagram is not reliable, as can be seen in the lower panel of Fig. 5. The cumulative plots do not reach plateaus around 1 fm, but instead tend to fall significantly up to 2 fm, or more. Once the cut moves beyond 1 fm it is no longer effective. The different behavior between the two stems from the different sampling strategies used for each [32]. Using the improved connected result, we find our final result for QED L , where the error is mostly statistical. We also include all systematic errors added in quadrature, including the hybrid O(a 2 ) error of the connected diagram. The systematic errors are summarized in Tab. III.

V. SUMMARY AND OUTLOOK
We have presented results for the hadronic light-by-light scattering contribution to the muon g − 2 from Lattice QCD+QED calculations with all errors under control. Large discretization and finite volume corrections are apparent but under control, and the value in the continuum and infinite volume limits is compatible with previous model and dispersive treatments, albeit with a large statistical error. Despite the large error, which results after a large cancellation between quark-connected and disconnected diagrams, our calculation suggests that lightby-light scattering can not be behind the approximately 3.7 standard deviation discrepancy between the Standard Model and the BNL experiment E821. Future calculations will reduce the error significantly. The calculations presented here strengthen the much anticipated test of the Standard Model from the new experiments at Fermilab and J-PARC, with the former planning to announce first results near the beginning of 2020.
Office of Science Facility supported under Contract DE-AC02-06CH11357. We gratefully acknowledge computer resources at the Oakforest-PACS supercomputer system at Tokyo University, partly through the HPCI System Research Project (hp180151, hp190137), the BNL SDCC computer clusters at Brookhaven National Lab as well as computing resources provided through USQCD at the Brookhaven and Jefferson Labs.

A. QED Test with different fitting forms
We test various fitting formulas from lattice calculation of the muonic leptonic light-by-light contribution to muon g − 2. The analytic result is known to be a µ = 0.371 × (α/π) 3 = 46.5 × 10 −10 using the conventional perturbative calculation.
In Ref. [32], we have performed the lattice calculation with three different lattice volumes, and each with three different lattice spacings. The results are listed in Tab. IV. Then, the continuum limit is calculated for each lattice volume and then extrapolate to infinite volume limit. In this paper, we adopt a different strategy. We use one formula which include both discretization effects and finite volume effects to fit all the data points. In the fitting forms listed below, we studied   • fit-product-form-3L-4a This is the form used in the QED light-by-light calculation in Section IV. The results of the fit are listed in the following The results of the fit is listed in the following The results of the fit is listed in the following table: Variable Value Statistical Error The results of the fit is listed in the following −c 1 (m µ a) 2 The results of the fit is listed in the following We discuss how we estimate all systematic errors presented in our results, including: strange quark contribution to the connected diagrams, and sub-leading disconnected diagrams' contribution. Many of the above systematic errors are resulting from lacking the corresponding terms in the fitting formula, Eq. (5), which we shall refer to as the "fit-plus-form-2L-2a-2ad-4ad". We therefore estimate the systematic errors by fitting with different fitting formulas with these terms included. We listed all the fitting forms below. • fit-product-form-2L-2a-2ad-4ad-lna Again, in these formulas, a I is the lattice spacing for the Iwasaki ensembles, 48I and 64I. We define it to be zero for I-DSDR ensembles. Similarly, a D is the lattice spacing for I-DSDR ensembles, and it is zero for Iwasaki ensembles. With this notation, the lattice spacings for all our ensembles are always equal to a = a I + a D . In two of the fitting forms, α S (a) is used to estimate the size of the a 2 log(a 2 ) contribution. We use the 2-loop α MS S with Λ QCD = 325 MeV: All these fit forms correct the O(a 2 ) lattice artifacts presented in DWF lattice calculations for Iwasaki ensembles (48I and 64I). The results obtained has a relative large statistical error mostly because of the continuum extrapolation. In Section IV, we introduced the hybrid continuum limit for the connected diagrams' contribution, where we calculate the continuum limit for contributions from the region r ≤ 1 fm and use the relatively precise 48I ensemble results without extrapolation for the r > 1 fm region. This hybrid procedure is possible due to the small size of the contribution from the region r > 1 fm. This is expected due to the exponential suppression from QCD mass-gap, and, in addition, for the connected diagrams, we label the three vertex locations on the quark loop that connect to the internal photons as x, y, z and require r = |x − y| ≤ min(|x − z|, |y − z|). This is a direct consequence of Eq. (3) in Ref. [33]. Operationally, this hybrid continuum limit for the connected diagrams is implemented by replacing the long distance region (r > 1 fm) of the 64I ensemble data with the corresponding 48I ensemble data: We perform the above replacement under the super jackknife procedure before the global fit together with all the results from other ensembles.
To estimate the hybrid O(a 2 ) systematic error resulting from the above replacement, we perform a slightly different replacement: The difference of the results obtained with these two replacements is used as the estimation of the hybrid O(a 2 ) systematic error.
We then estimate systematic errors due to lacking higher order finite volume or discretization terms in our default fitting form fit-plus-form-2L-2a-2ad-4ad by comparing the results between the fitting formulas: For the systematic errors due to not including (a) the strange quark contribution to the connected diagrams and (b) sub-leading disconnected diagrams' contributions, we performed lattice calculation with the 24D ensemble using the QED ∞ method [47]. The results are plotted in Fig. 8 and Fig. 9. We estimate the systematic error due to omitting these two contributions to be: ±0.3 × 10 −10 and ±0.5 × 10 −10 respectively. We then add all these systematic effects in quadrature as our final systematic error. We also use these procedure per jackknife sample and calculate the statistical uncertainty of the systematic error. This procedure is repeated for each of the "con", "discon", and "tot" contributions. The "tot" contribution is calculated by perform a single fit with the sum of "con" and the "discon" contributions for each ensemble. Since we only calculated the "con" contribution for the 48D ensemble, 48D is not included in the fit for "tot". As such, the "tot" fit is usually not exactly equal to the sum of the individual "con" and "discon" fits, which we label them as "sum". The results are shown in the Tab. V. Statistical errors for central values and systematic errors are listed in the table. We use the "tot" in this table as our central value.
We did the same exercise for fit-product-form-2L-2a-2ad-4ad as a comparison. The results are shown in Tab. VI. . Sub-leading disconnected diagrams' contribution. We only include the second diagram in Fig. 3. Only light quark contribution is calculated except for the tadpole part, where we reuse the calculation for the disconnected diagrams in HVP calculations described in Ref. [3] and both light quark and strange quark contribution is included. The remaining diagrams are equally or more suppressed.
Based on the experiences from the QED test, result from the fit-product-form-2L-2a-2ad-4ad fit the data better. However, in QCD calculations, fit-plus-form-2L-2a-2ad-4ad leads to smaller statistical error. At present level of accuracy, smaller statistical error is more important than potentially smaller systematic error. In particular, the systematic error is largely cancelled between the connected diagrams and the disconnected diagrams, while the statistical error does not. The fitting results with all the fitting forms are also summarized in Tab VII.
We have also performed the calculation without the hybrid continuum limit. The results are shown in Tabs. VIII, IX, X.
Finally, the detailed results and plots for each fitting forms with or without the hybrid continuum limit are listed in remaining tables and plots.