Hadronic light-by-light scattering in the anomalous magnetic moment of the muon

Hadronic light-by-light scattering in the anomalous magnetic moment of the muon $a_\mu$ is one of two hadronic effects limiting the precision of the Standard Model prediction for this precision observable, and hence the new-physics discovery potential of direct experimental determinations of $a_\mu$. In this contribution, we report on recent progress in the calculation of this effect achieved both via dispersive and lattice QCD methods.


Introduction
The magnetic moment of the muon is one of the most precisely measured quantities in particle physics. In units of e 2mµ · 2 , its value is given by the gyromagnetic factor g. The prediction that g = 2 was an early success of the Dirac equation, applied to the electron. The relative deviation of the gyromagnetic factor from the Dirac prediction is conventionally called anomalous magnetic moment, and is denoted by a µ ≡ (g − 2) µ /2. Remarkably, the quantity a µ has been directly measured to 0.54ppm of precision [1]. The Standard Model (SM) prediction for a µ , see e.g. [2], is currently at a similar precision level, 0.37ppm [3]. The precision of the SM prediction is entirely limited by the hadronic contributions. Specifically, the hadronic vacuum polarisation, which enters at O(α 2 ), and the hadronic light-by-light contribution a hlbl µ , which is of order α 3 , contribute in comparable amounts to the absolute uncertainty; their respective depiction as Feynman diagrams is shown in Fig. 1. The new E989 experiment at Fermilab is underway (see [4] and the presentation of A. Driutti at this conference), with the stated goal of improving the precision of the measurement by a factor four, and the E34 experiment at J-PARC (see [5] and the presentation of T. Mibe at this conference) plans to achieve a similar precision with a very different technique. It is therefore essential to improve the precision of the predictions for the hadronic contributions in order to enhance the new-physics sensitivity of the upcoming experimental results.
In view of the observations above, the theory precision requirements for the shortterm future are the following: for the hadronic vacuum polarisation contribution, a hvp µ ≈ 6900 · 10 −11 , the goal is to consolidate the currently quoted [6] precision of 0.35% obtained using a dispersive representation with experimental e + e − data as input, and to approach that level of precision in lattice calculations [7]; while for a hlbl µ ≈ 100 · 10 −11 a precision of 10 to 15% suffices. Clearly, both tasks are very challenging. In this talk, I focus on a hlbl µ and refer the reader to the talk of Ch. Lehner for the status of a hvp µ . The activities linked to the determination of a hlbl µ can be divided into four classes: 1. Model calculations, which constituted the only approach until 2014, are based on pole-and loop-contributions of hadron resonances, in some cases also on constituent quark loops.
2. Dispersive approaches allow one to identify and compute individual hadronic contributions in terms of physical observables, such as transition form factors and γ * γ * → ππ amplitudes.

3.
A dedicated experimental program is needed to provide input for the model & dispersive approaches, e.g. (π 0 , η, η ) → γγ * at virtualities Q 2 3 GeV 2 ; there is in particular an active program at BES-III on this theme, see the talk by Y. Guo at this conference.
4. In terms of lattice calculations, two groups (RBC/UKQCD and Mainz) have been working on formulating and carrying out a direct lattice calculation of a hlbl µ .
An important question is ultimately, how well the findings from the different approaches fit together. We begin by reviewing aspects of the model calculations and describing how one can test the assumptions underlying them; carry on to describe the status of the dispersive approaches and finally discuss in more detail several aspects of the lattice calculations of a hlbl µ .
As compared to earlier estimates, the pole contribution of the axial-vector mesons has been revised and is much smaller. Nevertheless, the central value of the estimate has changed little since the 2009 'Glasgow consensus' estimate of (105 ± 26) · 10 −11 [9]. Beyond the numerical result of the model calculations, it is worth recording some of the physics lessons learnt from them [9]: • A heavy (charm) quark loop makes a small contribution where Q c and m c are respectively the charm quark charge and mass, N c = 3 is the number of colors, m µ is the muon mass and α ≈ 1/137 is the fine-structure constant.
• For light-quarks, the most relevant degrees of freedom are the pions. The leading contribution in chiral perturbation theory, namely the charged-pion loop calculated in scalar QED, depends only on m µ /m π , with Numerically, this contribution amounts to a hlbl µ ≈ −45 · 10 −11 for the physical value of m µ /m π . Secondly, the neutral-pion exchange is positive and sensitive to the confinement scale [10,11], We note that, the pion decay constant F π ≈ 92 MeV being of order N 1/2 c , the contribution (3) is enhanced by a factor N c relative to the pion loop, Eq. (2). On the other hand, the latter is dominant in the limit m π → 0 with m µ /m π fixed. The ρ meson mass appears as the hadronic scale regulating an ultra-violet divergence, which appears if one assumes a virtuality-independent π 0 → γ * γ * coupling.
• For real-world quark masses, using form factors for the mesons is essential in obtaining quantitative results, and resonances up to 1.5 GeV can still be relevant. This makes a hlbl µ sensitive to QCD at intermediate energies, which is difficult to handle by analytic methods.
• Some information can be obtained from the operator-product expansion. Two closeby vector currents 'look like' an axial current from a distance. For that reason, the doubly-virtual transition form factors of 0 −+ and 1 ++ mesons only fall off like 1/Q 2 . This singles out the poles associated with pseudoscalar and axial-vector mesons as being particularly relevant. However, the coupling of an axial-vector meson to two real photons  [15].
is forbidden by the Yang-Landau theorem [12,13], suggesting that the pseudoscalar mesons π 0 , η, η are the most important pole contributions in a hlbl µ due to their unsuppressed coupling to two photons, in addition to their relatively light masses.
The applicability of the hadronic model to a hlbl µ can be tested by predicting the relevant four-point correlation function of the electromagnetic current j µ = f =u,d,s,... Q fq γ µ q and confronting the prediction with non-perturbative lattice QCD data. The Euclidean momentum-space four-point function at spacelike virtualities can indeed be computed in lattice QCD [14,15], and projected to one of the eight forward γ * γ * → γ * γ * scattering amplitudes, for instance In this particular case, the projectors R µν project onto the plane orthogonal to the vectors Q 1 and Q 2 and M TT thus corresponds to the amplitude involving transversely polarized photons. Dispersive sum rules have been derived for the forward amplitudes [16,17]. With ν = 1 2 (s+Q 2 1 +Q 2 2 ), a crossing-symmetric variable parametrizing the center-of-mass energy √ s, we can write a subtracted dispersion relation, where σ J corresponds to the total cross-section for the photon-photon fusion reaction γ * γ * → hadrons with total helicity J.
While experimental data exists for the fusion of real photons into hadrons, no such data is available for spacelike photons. In order to model the corresponding cross-section, we note that the contribution of a narrow meson resonance is It is then interesting to test whether all eight forward LbL amplitudes obtained from lattice computations can be described by such a sum of resonances via the dispersive sum rule. Essential ingredients in this parametrization of σ γ * γ * →hadrons are the transition form factors , describing the coupling of the resonance to two virtual photons. In the case of the neutral pion, a dedicated lattice QCD calculation of F π 0 γ * γ * was performed [18], thus allowing for a definite prediction for this contribution. For the other included hadronic resonances, which have quantum numbers J P C = 0 ±+ , 1 ++ , 2 ++ , a monopole or dipole parametrization of the virtuality-dependence of the transition form factors was chosen and fitted to the lattice data for the forward LbL amplitudes. In addition to the resonances, the Born expression for σ γ * γ * →ππ was included in the cross-section. A satisfactory description of the data was obtained in this way; see Fig. 2.
The calculation of the four-point correlation function of a quark bilinear in lattice QCD requires computing the Wick contractions of the quark fields, since the action is bilinear in these fields. The major difference with perturbation theory is that the quark propagators have to be computed in a non-perturbative gauge field background, which means inverting the sparse matrix of typical size 10 8 × 10 8 representing the discretized Dirac operator, on a source vector. The back-reaction of the quarks on the gauge field is taken into account in the importance sampling of the gauge fields. Five classes of Wick contractions contribute to the full four-point correlation function, as illustrated in Fig. 3. While the fully connected class of diagrams can be computed cost-effectively using 'sequential' propagators, the other classes require the use of stochastic methods. In [15], only the first two classes, denoted by the symbols (4) and (2+2), were computed, because the other three classes (3+1), (2+1+1) and (1+1+1+1) are expected to yield significantly smaller contributions. If this expectation is correct, and if the LbL amplitude is dominated by resonance exchanges, one can infer with what weight factors the isovector and the isoscalar resonances contribute to the leading contraction topologies (4) and (2+2). The isoscalar resonances contribute (with unit weight) to the class (2+2); the isovector resonances overcontribute with a weight factor 34/9 to class (4), while the (2+2) contractions compensate with a weight factor of −25/9 [19]. These counting rules have been used in describing the lattice data in Fig. 2. In particular, the large-N c inspired counting rules suggest that there is a large cancellation between the isovector resonances and the isoscalar resonances in the (2+2) class of diagrams, with the exception of the pseudoscalar mesons, due to the large mass difference between the π 0 and the η meson. Therefore, the contribution of the (2+2) diagrams to the light-by-light amplitudes was modelled as the η contribution, minus Thus the exploratory study [15] found that the LbL tensor (5) at moderate spacelike virtualities can be described by a set of resonance poles, much in the same way that a hlbl µ is obtained in the model calculations. It would be worth exploring this avenue further, in particular by increasing the precision of the lattice calculation.
3 Dispersive approach to a hlbl µ and its input Several dispersive approaches have been proposed to handle the complicated physics of hadronic light-by-light scattering [20][21][22]. Here we will mainly review aspects of the 'Bern approach' [21], which is the furthest developed at this point in time. It was shown that the full hadronic light-by-light tensor can be decomposed into 54 Lorentz structures [23], The Lorentz-invariant coefficients Π i are entirely determined by seven functions of the invariants q i ·q j combined with crossing symmetry. The 54 Lorentz structures are redundant, but they allow one to avoid kinematic singularities. The HLbL contribution to (g − 2) µ is then computed using the projection technique, i.e. directly at q = 0: TheΠ i are linear combinations of the Π i appearing in Eq. (9). Performing all "kinematic" integrals using the Gegenbauer-polynomial technique after performing a Wick rotation, the expression can be reduced to a three-dimensional integral, the master relation in this approach (τ = Q 1 · Q 2 /(|Q 1 | |Q 2 |)).
The contribution of the pole contributions associated with pseudoscalar mesons was worked out explicitly and clarified the way that the corresponding transition form factors are to be applied in this framework; see subsection 3.1 below. As a further result obtained as part of this approach, it was shown [24] that certain contributions in the dispersive approach to the pion loop could be handled rather accurately, The rescattering effects of the pions are being worked out for partial waves ≤ 2 [25]; first results by the Bern group for the s-wave were presented at the (g − 2) Theory Initiative Workshop [26]. An independent analysis of the γγ * → ππ process has also appeared very recently [27].

The transition form factor of the pion
The field-theoretic definition of the transition form factor of the pion involves a timeordered product of two vector currents, with p = q 1 + q 2 . A detailed dispersive analysis of the π 0 → γ * γ * transition form factor has recently been carried out [28], leading to the rather accurate result a HLbL,π 0 µ = 62.6 +3.0 −2.5 · 10 −11 .
In addition to experiments, important input for the dispersive approaches can be provided by lattice QCD. A first calculation of F πγ * γ * (q 2 1 , q 2 2 ) was carried out in lattice QCD with the two lightest quark flavors [18] and used to calculate a HLbL,π 0 µ , obtaining the result (65.0 ± 8.3) · 10 −11 . A model parametrization of the transition form factor was used here which incorporates known constraints at asymptotically large virtualities 2 . A second calculation in QCD, including also the dynamical effects of the strange quark, and using a more model-independent conformal-mapping parametrization of F πγ * γ * (q 2 1 , q 2 2 ), obtained the preliminary result a hlbl µ | π 0 = (60.4 ± 3.6) · 10 −11 [29]. The lattice and dispersive results are thus in excellent agreement, and comparable in precision. It is somewhat surprising how close the central values are to older estimates based on the simplest vector-meson dominance model of the form factor, e.g. a hlbl µ | π 0 ,VMD = 57.0 · 10 −11 [10]; however, the uncertainty of the result is now much better known.

The direct lattice calculation of HLbL in
The idea to directly calculate a hlbl µ in lattice QCD was pioneered in [30]. At first, the task was thought of as a combined QED + QCD calculation. Today's viewpoint is that the calculation amounts to a QCD four-point function, to be integrated over with a weighting kernel which represents all the QED parts, i.e. muon and photon propagators. Two collaborations have so far embarked on this challenging endeavour [7].
The RBC/UKQCD collaboration has performed calculations of a hlbl µ using a coordinatespace method in the muon rest-frame. The photon and muon propagators are either computed on the same L × L × L × T torus as the QCD fields -this approach goes under the name of QED L and was first published in [31]; or they are computed in infinite volume in a method called QED ∞ [32].
The Mainz group has used a manifestly Lorentz-covariant QED ∞ coordinate-space strategy, presented in [33], averaging over the muon momentum using the Gegenbauer polynomial technique. This technique relies on the anomalous magnetic moment of the muon being a Lorentz scalar quantity; a fact that has been used extensively in the phenomenology community.
A theoretical advantage of using the QED ∞ formulation is that no power-law finitevolume effects appear, which arise in QED L due to the massless photon propagators 3 . Specifically, in the (Euclidean) notation used by the Mainz group, the master equation for computing a hlbl µ is a hlbl The QED kernelL [ρ,σ];µνλ (x, y) is computed in the continuum and in infinite-volume; it consists of the muon and photon propagators depicted in Fig. 4. Once the two tensors appearing in Eq. (15) are contracted, the integrand is a Lorentz scalar, and the integrals over x and y reduce to an integral over three invariant variables, e.g. (x 2 , x · y, y 2 ). In this sense, Eq. (15) is the coordinate-space analogue of Eq. (11). In practical lattice calculations, one may carry out the four-dimensional summation over one of the vertices (say x) in full, because this tends to have a beneficial averaging effect. For fixed vertex position y, the four-dimensional summation over x can be arranged to be performed exactly on every gluon field configuration with an affordable number of operations for the most important Wick contraction classes (4) and (2+2). Sampling all values of the vertex position y would be computationally too costly; reducing instead the integral over y to a one-dimensional integral [33] allows one to sample the integrand reliably.
In order to get an idea of the length scales involved in the problem, one can inspect the coordinate-space dependence of the integrand for the pion pole [34]. An example is shown in the left panel of Fig. 5 for the physical pion mass. The curves correspond to different kernels, which are equivalent in infinite volume. They differ by x-or y-independent terms which do not contribute to the integral, but modify the integrand. In [32], where these subtraction terms were first introduced, it was shown that they can drastically reduce lattice discretization errors by forcing the kernel to vanish when some of the vertices coincide. From the left panel of Fig. 5, it is clear that the integrand is rather long-range in all three cases shown.

Status of lattice results
The calculation of hadronic light-by-light scattering in (g − 2) µ is still work in progress. The most recent refereed publication by the RBC/UKQCD collaboration [35] uses the QED L formulation and presents results for the Wick-contraction classes (4) and (2+2) on a 48 3 × 96 lattice at the physical pion mass with a lattice spacing of a −1 = 1.73 GeV. The contribution of the fully connected class of diagrams is a HlbL (4) µ = (116.0 ± 9.6) × 10 −11 , while the (2+2) diagrams yields a HlbL (2+2) µ = (−62.5±8.0)×10 −11 . Together, they amount to [35] a hlbl µ = (53.5 ± 13.5) · 10 −11 , where the quoted error is statistical. This represents the first lattice result for the two leading Wick-contraction topologies. However, the authors acknowledge that "The finitevolume and finite lattice-spacing errors could be quite large and are the subject of ongoing research." We note that the total is about a factor two lower than the model estimates. A possible explanation is that in the latter, neglected contributions could be more important  (2) , which vanishes by construction whenever x or y vanishes. The integrals over the vertices in x and z (see Fig. 4) have already been performed at this stage. For comparison, the corresponding (continuum, infinite-volume) integrand for the π 0 pole contribution to a hlbl µ is displayed, multiplied by an enhancement factor of 3 (lower edge of the band) and 34/9 (upper edge of the band) to account for the absence of the disconnected diagrams; see the end of section 2 for an explanation of these weight factors. than so far expected; or, since the QED L method has O(1/L 2 ) finite-size effects, the latter could be responsible for a large systematic error. It was noted [15] that, based on the model estimate and large-N c inspired arguments, one would expect a HlbL (2+2) µ to be dominated by the (π 0 , η, η ) exchange and approximately −150 · 10 −11 in size; the fully connected diagrams would then have to amount to about 250 · 10 −11 in order to recover the model result a hlbl µ ≈ 100 · 10 −11 discussed in section 2. L. Jin presented an update of the RBC/UKQCD calculation at this year's lattice conference. An extrapolation to infinite volume and zero lattice spacing based on several ensembles yielded the results a HlbL (4) µ = (282 ± 40) · 10 −11 and a HlbL (2+2) µ = (−163 ± 34) · 10 −11 , resulting in the sum a hlbl µ = (119 ± 53) · 10 −11 . These values are much more in line with the expectation from the model calculations, but the extrapolation, and the cancellation between the two Wick-contraction topologies, enhances the relative error on the final result.
Both the RBC/UKQCD collaboration and the Mainz group have started generating lattice results with their respective QED ∞ method. As a very recent development, the RBC/UKQCD collaboration has performed a calculation of the three leading diagram topologies (4), (2+2) and (3+1) on a coarse lattice at the physical pion mass; see Fig. 6. The (3+1) topology is found to make a negligible contribution [36]. The Mainz group has computed and analyzed the fully connected set of diagrams (4) [37]. It uses rather fine lattices, with a typical lattice spacing of a = 0.064 fm, on the other hand the simulated pion masses (m π 200 MeV) lie above the physical value. The integrand obtained at m π = 340 MeV is displayed in the right panel of Fig. 5, and compared to the prediction corresponding to the neutral pion pole, including the appropriate enhancement weight factor, as discussed at the end of section 2. The pion-pole integrand, predicted with a VMD transition form factor, provides a surprisingly good approximation to the lattice data. as a function of its upper limit, the integration variable being the maximum distance between any two internal vertices. The QED ∞ formulation including subtractions [32] was used. Separately (from left to right) for the Wick-contraction classes (4), (2+2) and (3+1). In the latter case, only those Wick-contractions are included in which the external photon is attached to the loop containing three vertices.
Tests of the systematic errors on the fully connected set of diagrams at m π = 280 MeV are shown in Fig. 7: finite-size effects and discretization effects appear to be under control. It remains to be seen how well controlled the extrapolation in the pion mass will be. Restricting the sums over the vertex positions in a systematically controlled way to the regions that contribute appreciably is likely to reduce the statistical noise.

Conclusion
The hadronic light-by-light scattering contribution is, along with the hadronic vacuum polarization contribution, the leading source of uncertainty in the Standard Model prediction for the anomalous magnetic moment of the muon, (g − 2) µ . For decades, a framework which offers a systematically improvable prediction was lacking. This has now changed, with significant progress having been made in the dispersive as well as lattice QCD approaches. Even though model calculations will soon become superseded, valuable lessons have been learnt from them, which can help control the systematic errors in the ab initio approaches.
Given that the quantity a hlbl µ involves three spacelike and one quasi-real photon, lattice QCD is in a good position to provide a first-principles prediction -no analytic continuation is required. At the same time, dealing with the four-point function of the vector current is pushing the field into a territory on which the community has little prior experience, hence dealing with the statistical and the systematic errors, such as the finite lattice spacing and finite-volume errors, requires special attention. A cross-check between at least two independent lattice collaborations is hence extremely valuable, as is a comparison of the lattice results with the results of the dispersive approaches. In the latter case, it will be especially interesting to see how the dispersive treatment of ππ intermediate states differs quantitatively from previous estimates based on a narrow (scalar and tensor) resonance exchange approximation. as a function of the upper limit of the integration variable |y| at a pion mass of m π = 280 MeV. The integrals over the vertices in x and z have already been performed at this stage. The kernel used is L (2) , which vanishes by construction whenever x or y vanishes. The left panel compares results from two different lattice volumes; the right panel compares two different lattice spacings.