Lattice calculation of the hadronic leading order contribution to the muon $g-2$

The persistent discrepancy of about 3.5 standard deviations between the experimental measurement and the Standard Model prediction for the muon anomalous magnetic moment, $a_\mu$, is one of the most promising hints for the possible existence of new physics. Here we report on our lattice QCD calculation of the hadronic vacuum polarisation contribution $a_\mu^{\rm hvp}$, based on gauge ensembles with $N_f=2+1$ flavours of O($a$) improved Wilson quarks. We address the conceptual and numerical challenges that one encounters along the way to a sub-percent determination of the hadronic vacuum polarisation contribution. The current status of lattice calculations of $a_\mu^{\rm hvp}$ is presented by performing a detailed comparison with the results from other groups.


Introduction
The anomalous magnetic moment of the muon, a µ , has been measured by the E821 experiment at BNL [1] with a total precision of 0.54 ppm, i.e. a µ ≡ 1 2 (g − 2) µ = 116 592 080(54)(33) · 10 −11 . (1) While the theoretical prediction based on the Standard Model (SM) is known at the same level of precision, it differs by 3.5 standard deviations [2]. This discrepancy is currently one of the most intriguing hints for the possible existence of physics beyond the SM. While the absolute value of the theoretical estimate of a µ is dominated by the contributions from QED, with only a small fraction coming from strong and weak interaction effects, the situation is completely reversed regarding the assigned uncertainty: Indeed, the dominant errors arise from the strong interaction, in particular from the hadronic vacuum polarisation (HVP) and hadronic light-by-light scattering (HLbL) contributions, a hvp µ and a hlbl µ , respectively. New experiments, which are designed to reduce the error even further, are either being prepared (E34 at J-PARC [3]) or have started data taking already (E989 at Fermilab [4]). Since the design sensitivity of both experiments exceeds the theoretical uncertainty by far, it is of utmost importance to achieve a significant reduction of the latter in the coming years.
Simulations of QCD on a space-time lattice have emerged as a versatile tool for the treatment of the strong interaction in the non-perturbative regime from first principles. Lattice QCD provides accurate estimates for SM parameters (strong coupling α s , quark masses), decay constants and form factors relevant for flavour physics, quantities describing structural properties of the nucleon and many other observables of phenomenological relevance [21]. Recent years have also seen enormous progress regarding calculations of the hadronic contributions to the muon g − 2 (see Ref. [22] for a review).
A first-principles approach to the muon g − 2 based on lattice QCD avoids both the reliance on experimental data, as well as model estimates. However, in order to have an impact on tests of the SM, lattice QCD must deliver results for a hvp µ with a total error well below 1% and a determination of a hlbl µ at the 10 − 15% level. Whether or not this can be achieved is also the subject of a planned White Paper summarising the current status of the theoretical prediction ahead of the first results from E989.

Lattice formalism and general issues
The hadronic vacuum polarisation contribution can be expressed as a convolution integral involving the subtracted vacuum polarisation function, i.e. [23,24] where the integration is over Euclidean momenta, Π(Q 2 ) is obtained from the vacuum polarisation tensor and J µ = 2 3ū γ µ u− 1 3d γ µ d− 1 3s γ µ s+. . . is the electromagnetic current. The evaluation of the convolution integral using lattice QCD data for the correlator is hampered by the fact that the weight function f (Q 2 ) is strongly peaked near the muon mass, i.e. for |Q| ∼ m µ . Since |Q| is inversely proportional to the spatial lattice extent L, the peak region corresponds to lattice volumes that are much larger than those that can be realised in current simulations. Thus, there is a lack of lattice data in the region where the integral receives its dominant contribution.
An alternative formulation in position space, suggested by Bernecker and Meyer [25], is known as the timemomentum representation (TMR) where G(t) denotes the spatially summed vector-vector correlator. Empirically, one finds that the region t 3 fm contributes ≈ 3% to the integral in Eq. (4). In order to determine a hvp µ with a precision of 1% or better, it is clear that G(t) must be calculated with sufficient accuracy in the long-distance regime. This, however, is a challenging task, since the noise-to-signal ratio in G(t) increases exponentially with Euclidean time t.
While both integral representations are equivalent in the sense that they are formulated in terms of the currentcurrent correlator, the TMR turns out to be technically simpler: this is due to the fact that detailed spectral information is available for the vector correlator in the longdistance regime which is dominated by the isovector channel. Furthermore, the calculation of quark-disconnected diagrams, which requires major computational resources (see Section 3.5 below) is technically more straightforward.
The main challenges facing lattice QCD regarding a calculation of a hvp µ with a precision that will have a decisive impact on SM tests are the following: First, the statistical error of the evaluation of Eq. (4) should lie significantly below 1%. This implies that the vector correlator G(t) must be accurately constrained for t 3 fm. The long-distance regime of G(t) receives significant finitevolume corrections for typical lattice sizes, i.e. for m π L 4, which have to be quantified. The isoscalar part of G(t) requires accurate knowledge of the contributions from quark-disconnected diagrams which typically exhibit a high level of statistical noise. Finally, one must include the effects of strong (m u m d ) and electromagnetic isospin breaking corrections. In the next section we will describe how these issues are addressed in our recent calculation [26][27][28].
3 Lattice calculation of a hvp µ Our calculations employ the gauge ensembles generated by the CLS consortium, using N f = 2 + 1 flavours of O(a) improved Wilson quarks, at four different values of the lattice spacing, a range of pion masses down to the physical one and volumes satisfying m π L ≥ 4 throughout. For a description of technical details, including the definition of the vector currents used to compute the correlator G(t), we refer the reader to Refs. [26,28,29].

Controlling the infrared regime
The correlation function G(t) of the electromagnetic current defined in Eq. (4), can be split into a sum over the quark-connected contributions from individual flavours plus the contribution from quark-disconnected diagrams This expression holds in the isosymmetric case, m u = m d ≡ m l , with the numerical prefactors arising from the electric charges in the definition of the electromagnetic current. For the following discussion it is convenient to consider the isospin decomposition of the vector correlator, i.e.
where G I=1 (t) and G I=0 (t) denote the correlators in the isovector and isoscalar channels (see Section 2.2.3 in [22] for further details). Since the spectral function in the isoscalar channel vanishes below the 3-pion threshold, one expects G(t) to be dominated by the lowest-energy state in G I=1 (t) as t → ∞. We have studied three different and complementary methods to accurately constrain the longdistance regime of the isovector correlator G I=1 (t): The first makes use of a dedicated calculation of the energy levels in the isovector channel, by observing that where the energies ω n are related to the scattering momentum k and the pion mass via ω n = m 2 π + k 2 [30,31]. Figure 1 shows that the 3-4 lowest-lying states saturate the correlator for distances t 2 fm [26].
The second procedure, known as the "bounding method" [32,33], exploits the positivity of the spectral sum for G I=1 (t), which implies that  bounds on G(t) can be significantly sharpened using the knowledge from the dedicated calculation of G I=1 (t) [26]. By subtracting the N − 1 lowest-lying states one defines the correlatorG(t) viã which satisfies 0 ≤G(t) ≤G(t cut ) e −ω N (t−t cut ) , where ω N is the energy of the Nth state. In Fig. 2 we show a comparison of the ordinary and improved bounding methods applied to the determination of the (quark-connected) light quark contribution to a hvp µ , as a function of t cut . One finds thatG(t) saturates the bounds implied by the improved bounding method already for t cut 1.5 fm, which, when averaged over a few timeslices, leads to a stable and statistically precise estimate of a hvp µ . The third method to constrain the long-distance regime of the isovector correlator makes use of the information provided by the timelike pion form factor F π (ω) [34]. It provides a relation between F π (ω) and the amplitudes A n in Eq. (7), which requires knowledge of the Lellouch-Lüscher factor [35], which, in turn, is determined from the scattering phase shifts [30,31] (see Eq. (13) below). A detailed study of this method in two-flavour QCD can be found in Ref. [36].

Finite-volume effects
The long-distance regime of the vector correlator is closely related to the issue of finite-volume corrections. The latter can be determined by inserting the difference of the vector correlator in infinite and finite volume, i.e.
where G(t) is the correlator in finite volume, and G(t, ∞) is its infinite-volume counterpart. Realising that the longdistance regime is most relevant for the estimation of the finite-volume correction, one can focus on the isovector correlator. In infinite volume, its spectral representation reads where the spectral function is given in terms of the pion form factor F π (ω) as In finite volume, the expression for the isovector correlator G I=1 (t) has been listed Eq. (7), with amplitudes given by where k is the scattering momentum. Knowledge of the pion form factor F π (ω) allows one to determine the differ- and thus the correction ∆a hvp µ . In our evaluation of finite-size effects, we have employed the Gounaris-Sakurai (GS) representation [37] for F π (ω), using the energy levels and amplitudes, (ω n , A n ), as input. Alternatively, we have used m ρ and the coupling g ρππ , taken from a dedicated calculation of F π (ω) on a subset of our ensembles [38] . We stress that the somewhat simplistic GS model is used only for the determination of the finite-size correction but not for the treatment of the long-distance tail of the correlator itself. In order to test the ability of this approach to provide an accurate description of finite-size corrections, we have compared the results for the integrand computed at a fixed pion mass of 280 MeV for two different volumes, corresponding to m π L = 3.8 and 5.8, respectively. Applying the correction to the data obtained on the smaller volume produces values that are in excellent agreement with those observed for m π L = 5.8. We conclude that finite-size effects are well described by the GS parameterisation of F π (ω).

Scale setting uncertainty
At first it may seem surprising that a dimensionless quantity such as a hvp µ should be affected by the lattice scale. However, the kernelK(t) in Eq. (4) is actually a function of (tm µ ) 2 , and hence it is necessary to express the argument t/a of the vector correlator in units of the muon mass. This requires the calibration of the lattice spacing a. Furthermore, the quark masses enter the determination of a hvp µ indirectly via their feedback onto the gauge field during the generation of the ensembles. Therefore, the HVP contribution a The relative error on the HVP contribution, ∆a hvp µ /a hvp µ , which arises from the scale setting uncertainty ∆Λ/Λ is given by where the ellipses represent derivatives with respect to M u , M d , . . .. The logarithmic derivative with respect to M µ can be worked out in terms of an integral representation: where J(t) = tK (t) −K(t). Using input for the correlator and the representation of the kernel function described in Appendix B of [39], the proportionality factor between ∆a hvp µ /a hvp µ and ∆Λ/Λ can be estimated to be ≈ 1.8. Hence, the scale setting uncertainty must, at least, be smaller by a factor two compared to the target precision in the determination of a hvp µ . The contributions arising from the quark mass dependence have also been quantified and turn out to be subdominant [39,40]. In our calculation, we use the lattice scale determined in [41], i.e. Λ −1 = √ 8t 0 = 0.415(4)(2) fm.

Isospin breaking
The effects from strong and electromagnetic isospin breaking must be taken into account to obtain a lattice QCD estimate for a hvp µ whose precision can compete with that of the data-driven dispersive approach. One way to include isospin-breaking effects is the so-called "stochastic method" which is based on the direct Monte Carlo evaluation of the path integral in the presence of electromagnetism. The effects of strong isospin breaking are incorporated either by choosing different quark masses, m u m d or via reweighting techniques. The second method to quantify the effects of isospin breaking is based on the perturbative expansion of the path integral in powers of the fine structure constant α = e 2 /4π and the mass splitting m d − m u [42,43]. A detailed comparison of both methods was performed in [44].
The Mainz group has started to determine isospinbreaking corrections to a hvp µ using the perturbative approach [27,45,46]. In this way, isospin-breaking effects can be extracted from additional correlators of the vector  current that involve insertions of the scalar density and explicit photon propagators (in Coulomb gauge), computed on our existing iso-symmetric gauge ensembles. Since the photon is a massless unconfined particle, one expects strong finite-volume effects in the presence of electromagnetism. Our calculations are based on the QED L prescription [47] in which the spatial zero modes of the photon field are set to zero. The necessary adaptation of QED L in the presence of open boundary conditions, which are employed in the generation of most of our gauge ensembles, is discussed in [46]. A list of the additional diagrams for the vector correlator, as well as first results on the renormalisation factors of the electromagnetic current and the renormalised hadronic vacuum polarisation function are shown in [27] and Fig. 3. Figure 4 shows a compilation of recent results for isospin-breaking corrections obtained by different collaborations. One finds that the typical size of isospin-breaking effects is (5 − 10) · 10 −10 , which corresponds to a onepercent correction to a hvp µ . While the separation of strong and electromagnetic isospin-breaking effects is, in principle, scheme-dependent, the results indicate that the latter are subdominant.  Figure 4. Compilation of recent results for the isospin-breaking correction δa IB µ to the hadronic vacuum polarisation, determined by ETMC [48], RBC/UKQCD [49], BMW [50] and Fermilab/HPQCD/MILC [51]. Solid circles denote the full (strong and electromagnetic) isospin-breaking correction, while open symbols represent strong isospin breaking only. The result marked by the green cross is based on a phenomenological estimate.

Quark-disconnected diagrams
While the correlator G(t) contains both quark-connected and quark-disconnected contributions, the latter are difficult to quantify with the desired level of precision, owing to the large computational cost when employing the conventional source method (see, for instance, the discussion in Section 2.3 of Ref. [22]). Therefore, stochastic methods are applied which, despite being numerically quite efficient, introduce stochastic noise in addition to the statistical fluctuations due to the Monte Carlo integration over the gauge ensemble. Several methods have been devised to combat stochastic noise, such as performing a hopping parameter expansion of the propagator [52,53], low-mode deflation [54], hierarchical probing and Hadamard vectors [55], as well as the use of a Lorentz covariant coordinate space formulation [56] instead of the spatially summed vector correlator.
An important observation, reported initially in [57], is the partial cancellation of the stochastic noise in the difference between the contributions of the light and strange quarks, l − s. This has been exploited in a recently proposed method [58] which uses the so-called one-end-trick [59,60] in the calculation of the l − s contribution, combined with the hopping parameter expansion applied to a single heavy quark flavour. Our initial runs show that this method outperforms hierarchical probing by more than an order of magnitude in the case of the vector-vector correlator.
We have computed the quark-disconnected contribution to G disc (t) (see Eq. (5)) on a subset of our gauge ensembles. During our calculation we have found that it is difficult to constrain the long-distance contribution of G disc (t) to the integrand, since the signal is lost for t 1.5 fm. We have therefore adopted another strategy, which applies the bounding method to the isoscalar correlator G I=0 (t) which contains the disconnected contribution. Neglecting the contributions from the charm quark, the positivity of the isoscalar spectral function implies that the isoscalar correlator satisfies [26]     where we have assumed m ω ≈ m ρ , since we do not have a dedicated calculation of the mass of the ω-meson at our disposal. The quark-disconnected contribution (a hvp µ ) disc to the HVP is then obtained by subtracting the connected contributions of the light (l) and strange (s) quarks from the isoscalar part, i.e.
In Fig. 5 we show the results for (a hvp µ ) disc as a function of (m 2 K − m 2 π ) 2 . After extrapolating the results to the physical mass, we obtain (a hvp µ ) disc = (−23.2 ± 2.2 ± 4.5) · 10 −10 , where the first error is statistical, and the second is estimated from the ambiguity in the ansatz for the chiral extrapolation. A compilation of recent results for the quarkdisconnected contribution is shown in Fig. 6.

Results at the physical point
Our results for the hadronic vacuum polarisation contribution have been obtained on ensembles spanning a range of lattice spacings and pion masses, including the physical mass. We have subjected our results for a . While the estimate of a hvp µ does not contain the correction due to isospin breaking, the size of the correction from Ref. [48] has been added in quadrature to the error.
to the physical point. To this end we have investigated several ansätze, taking guidance from Chiral Perturbation Theory [26]. In particular, for the light quark contribution we have used and compared ansätze that reproduce the expected singular behaviour as m π → 0.
The fact that we have employed two different discretisations for the vector current [61] allowed us to extrapolate the results simultaneously to a common value at the physical point. Only for the charm quark contribution did we observe large lattice artefacts for the data extracted from the local-local vector correlator, so that we decided to extrapolate the results from the local-conserved correlator, which show much milder discretisation effects. At the physical point we obtain (in units of 10 −10 ): where the first error is statistical and the second gives an estimate of the systematic uncertainty. The latter is dominated by the form of the ansatz for the chiral extrapolation and the scale setting uncertainty. Our final result is obtained by adding the quark-disconnected contribution from Eq. (19) to the sum of the above estimates. We obtain a hvp µ = (720.0 ± 12.4 stat ± 9.9 syst ) · 10 −10 .
Since our calculation of isospin-breaking effects has not been completed, we have increased the systematic error by adding in quadrature the result for the full isospin-breaking correction from Ref. [48]. The pie charts in Fig. 7 show the individual contributions to the estimate and the variance in Eq. (23). The total uncertainty of our result is 2.2%.

Comparison, conclusions and outlook
Our result [26], shown in Eq. (23), is compared to other recent estimates from lattice QCD [39,[48][49][50][62][63][64] and to the data-driven approach in Fig. 8.  [62]. Given that a range of different discretisations are employed and that the strange and charm quarks make only small contributions to the tail of G(t), the observed consistency can be taken as an indication that lattice artefacts are under good control.
By contrast, the results for the connected light-quark contribution, (a hvp µ ) l con , show a much higher degree of scatter even though there is a certain degree of correlation among different estimates: For instance, the calculations by Aubin et al. [64] and Fermilab-HPQCD-MILC (FHM 19) [63] share a subset of the same gauge ensembles. Furthermore, most groups employ a similar formalism to quantify finite-volume effects. A weighted average of the results for (a hvp µ ) l con yields χ 2 /d.o.f. = 14.62/6 2.44 which illustrates that the scatter among the results is large compared to the quoted total uncertainty. By contrast, weighted averages of the estimates for (a hvp µ ) c and (a hvp µ ) s yield χ 2 /d.o.f. ≈ 1 in either case. These observations call for further scrutiny of the light connected contribution, in particular regarding the estimation of the tail of G(t) as well as finite-volume effects.
Similar comparisons of isospin-breaking corrections and quark-disconnected contributions are shown in Figs. 4 and 6. While the agreement of isospin-breaking effects among different groups is quite good, the same cannot be said for the quark-disconnected contribution which, in the case of our own calculation, turns out to be a factor two larger in magnitude than what has been reported by other collaborations. Using a new technique that takes inspiration from Refs. [58][59][60], we are currently computing the quark-disconnected contribution with much higher statistical accuracy and over a wider range of ensembles than reported previously.
The comparison between lattice estimates and the result from dispersion relations, which is represented by the red vertical band, 1 shows that the overall precision of lattice calculations is not yet sufficient to distinguish between the data-driven approach and the "no new physics" scenario, which is represented by the green band in the right- (a hvp µ ) c · 10 10 Figure 8. Recent results for the hadronic vacuum polarisation contribution from Aubin et al. [64], Mainz/CLS [26], Fermilab/HPQCD/MILC (FHM) [63,65,66], PACS [62], ETMC [48,67,68], RBC/UKQCD [49] and BMW [50]. From left to right the panel represent the charm, strange, light (connected) and total contributions to a hvp µ . The vertical red line indicates the result from the data-driven analysis of Ref. [7]. The green vertical band corresponds to the "no new physics" scenario. most panel of Fig. 8. The latter corresponds to the difference between the experimental result for a µ and the SM prediction without the leading-order HVP contribution.
Further efforts are underway to reduce the overall uncertainty of lattice estimates for both the hadronic vacuum polarisation and hadronic light-by-light scattering contributions. The status, prospects and opportunities for lattice calculations applied to the muon (g − 2) will also be discussed in a White Paper which is currently in preparation.