Hadronic light-by-light scattering contribution to the muon g-2 on the lattice

We briefly review several activities at Mainz related to hadronic light-by-light scattering (HLbL) using lattice QCD. First we present a position-space approach to the HLbL contribution in the muon g-2, where we focus on exploratory studies of the pion-pole contribution in a simple model and the lepton loop in QED in the continuum and in infinite volume. The second part describes a lattice calculation of the double-virtual pion transition form factor F_{pi^0 gamma^* gamma^*}(q_1^2, q_2^2) in the spacelike region with photon virtualities up to 1.5 GeV^2 which paves the way for a lattice calculation of the pion-pole contribution to HLbL. The third topic involves HLbL forward scattering amplitudes calculated in lattice QCD which can be described, using dispersion relations (HLbL sum rules), by gamma^* gamma^* ->hadrons fusion cross sections and then compared with phenomenological models.


Introduction
The anomalous magnetic moment of the muon has served for many years as a precision test of the Standard Model [1][2][3][4] and it has also played an important role in many presentations at this meeting. There is a discrepancy of 3 − 4 standard deviations between experiment and theory for some time now. This could be a signal of New Physics [1,5], but the uncertainties in the theory prediction from hadronic vacuum polarization (HVP) and hadronic light-by-light scattering (HLbL) make it difficult to draw firm conclusions. In view of upcoming four-fold more precise new experiments at Fermilab and J-PARC [6], these hadronic contributions need to be better controlled.
The improvement for HVP looks straightforward, with more precise experimental data from various experiments on hadronic cross-sections as input for a dispersion relation [7]. But also lattice QCD is getting more and more precise [8], and, hopefully, also a new method using muonelectron scattering to measure the running of α and the HVP in the spacelike region [9], will be feasible at some point with the required precision.
On the other hand, the HLbL contribution to the muon g − 2, see Fig. 1, has only been calculated using models so far [1,4,10] and the frequently used estimates from Refs. [1,11,12] (revised slightly in Ref. [2]) suffer from uncontrollable uncertainties. In view of this, dispersion relations have been proposed a few years ago [13][14][15] (see also the very recent new proposal in Ref. [16]) to determine the presumably numerically dominant contributions from a single neutral pion-pole (light pseudoscalarpole) and from the two-pion intermediate state (pion-loop) based on input from experimental data for the dispersion * Speaker, e-mail: nyffeler@uni-mainz.de relations. Still some modelling will be needed to estimate the contributions from multi-pion intermediate states, like the axial-vector contribution.
Finally, lattice QCD was proposed some time ago as a model-independent, first principle approach to the HLbL contribution in the muon g − 2 and some promising progress has been achieved recently by the RBC-UKQCD collaboration [17].
Independently, also the lattice group at Mainz has studied HLbL in recent years. We used complementary approaches to tackle the full HLbL contribution in the muon g − 2 with a new position-space approach [18][19][20], studied the pion transition form factor (TFF) with two virtual photons on the lattice to evaluate the pion-pole contribution to HLbL [21] and we analyzed HLbL forward scattering amplitudes that can be compared using dispersion relations (HLbL sum rules) [22] to phenomenological models for photon fusion processes [18,23,24]. These three topics will be discussed in the following three sections. For all the details, we refer to the quoted papers. arXiv:1801.04238v1 [hep-lat] 12 Jan 2018 2 Position-space approach to HLbL in the muon g − 2 on the lattice The HLbL scattering contribution to the anomalous magnetic moment of the muon in Fig. 1 can be split into a perturbative QED kernelL that describes the muon and photon propagators and a non-perturbative QCD four-point function i Π (denoted by a blob in the Feynman diagram) that will be evaluated on the lattice. The projection on the muon g − 2 yields the master formula (in Euclidean space notation) [18,19] a HLbL with the spatial moment of the four-point function We evaluate the QED kernel in the continuum and in infinite volume and thereby avoid 1/L 2 finite-volume effects from the massless photons. Since Lorentz covariance is manifest in our approach, the eight-dimensional integral in Eq. (1) can be reduced, after contracting all indices, to a three-dimensional integral over the Lorentz invariants x 2 , y 2 and x · y.
The QED kernelL [ρ,σ];µνλ (x, y) can be decomposed into several tensors The G A δ[ρσ]µανβλ are traces of gamma matrices and just yield sums of products of Kronecker deltas. The tensors T (A) αβδ are decomposed into a scalar S , vector V and tensor T part These can be parametrized by six weight functions that depend on the three variables x 2 , x · y = |x||y| cos β and y 2 . The semi-analytical expressions for the weight functions have been precomputed to about 5 digits precision and stored on a three-dimensional grid. For illustration, the expression for the weight functionsḡ (2) has been given in Ref. [19]. In Fig. 2 we show the two weight functions g (1) andḡ (2) as a function of |x| < 12 fm, for a fixed value of |y| = 0.506 fm and three values of cos β. Plots for all six weight functions can be found in Ref. [20]. To test our semi-analytical expressions and the software for the QED kernel, we have computed the π 0 -pole contribution to HLbL in a simple vector-meson dominance (VMD) model as well as the lepton-loop contribution to a LbL µ in QED, where the results are well known. To this aim, we first derived analytical expressions for these contributions to the four-point function i Π ρ;µνλσ (x, y) in position-space, see Ref. [20] for details.
In figure 3 we plot the integrand f (|y|) of the final integration over |y| of the HLbL contribution a HLbL From the plot of the integrand for the pion-pole contribution in Fig. 3 one observes that this contribution to a HLbL µ is remarkably long-range with a long negative tail at large |y|. One expects an exponential decay ∼ e −cm π |y| of the correlation function. But this seems to be countered by some non-negligible power-like behavior |y| n . For pion masses m π = 300 − 900 MeV we reproduce the known results, which can be easily obtained from the three-dimensional integral representation in momentum space given in Ref. [1], at the percent level. On the other hand, for the physical pion mass, one will need rather large lattices of the order of 5−10 fm to capture the negative tail at large |y| in a QCD lattice simulation. Hopefully we can correct for finite-size effects on this contribution, by computing the relevant neutral pion transition form factor on the same lattice ensembles, see Ref. [21] and Section 3.
For the lepton-loop in QED the behavior of the integrand for small |y| is compatible with f (|y|) ∝ |y| log 2 (|y|). This is quite steep and means that we probe the QED kernel at small distances. In addition, the height of the positive peak grows with smaller masses m l of the lepton in the loop. Furthermore there is again a long negative tail at large |y| which demands the use of a large size for the grid where the weight functions have been calculated, in particular for a lepton in the loop that is lighter than the muon. For m l = m µ , 2m µ we reproduce the analytically known results for a LbL µ in QED [25] at the percent level, see Table 1. On the other hand, for the lightest lepton mass m l = m µ /2, some further refinements of our numerical evaluation are needed. Once this is achieved, we plan to make the QED kernel publicly available. Table 1. Results (×10 11 ), precision and deviation for the lepton-loop contribution to LbL in QED with our approach compared to the known results [25]. The first uncertainty originates from the three-dimensional numerical integration, the second from the extrapolation of the integrand to small |y|.   (1) (left) andḡ (2) (right) for |y| = 0.506 fm and three values of cos β. Apart from the scale, the weight functionsḡ (0) ,l (1) andl (2) have the same shape asḡ (1) , whereasl (2) looks similar toḡ (2) .

Lattice calculation of the pion transition
The pion transition form factor can be defined from the following correlation function in Euclidean space, using the methods proposed and used before in Ref. [26] provided the photon virtualities satisfy q 2 1,2 < min(M 2 ρ , 4m 2 π ) to avoid poles in the analytical continuation from the original definition of the form factor in Minkowski space to Euclidean space. The free real parameter ω 1 denotes the zeroth (energy) component of the four-momentum q 1 = (ω 1 , q 1 ).
The main object to compute on the lattice is the threepoint function Here τ = t i − t f is the time separation between the two vector currents and t π = min(t f − t 0 , t i − t 0 ). The matrix element in Eq. (10) with an on-shell pion is obtained by considering the limit of large t π . With the definitions one obtains The overlap factor Z π and the pion energy E π can be obtained from the asymptotic behavior of the pseudoscalar two-point function.
We choose the pion rest frame p = 0, with photons back-to-back spatially q 2 = − q 1 , where the kinematical range accessible on the lattice is given by q 2 . we obtain mostly spacelike photon virtualities up to |q 2 1,2 | ≈ 1.5 GeV 2 , as can be seen in Fig. 4 (left). In practice, discrete values of ω 1 have been used to sample the momenta. Note that on the lattice it is actually easier to access the double-virtual TFF F π 0 γ * γ * (q 2 1 , q 2 2 ), in particular near the diagonal q 2 1 = q 2 2 , than the single-virtual form factor F π 0 γ * γ * (q 2 1,2 , 0) along the two axis, in contrast to the situation in experiments [28].
From the theoretical side, the form factor is constrained by the chiral anomaly such that F π 0 γ * γ * (0, 0) = 1/(4π 2 F π ) [29] (in the chiral limit). For the single-virtual form factor one expects the Brodsky- [30]. The precise value of the prefactor is, however, under debate. The double-virtual form factor, where both momenta become simultaneously large, has been computed using the OPE at short distances. In the chiral limit the result reads In order to get a result for the double-virtual TFF F π 0 γ * γ * (q 2 1 , q 2 2 ) in the continuum and for the physical pion mass, we fit our lattice data, obtained for the eight ensembles with different lattice spacings a and pion masses m π , with three simple models: vector-meson dominance (VMD), lowest-meson dominance (LMD) and LMD+V, see Refs. [32,33] for details about these models. The often used VMD model fulfills the BL behavior for the singlevirtual case (and describes quite well the available experimental data below 2 − 3 GeV 2 ), but falls off too fast for the double virtual case F VMD π 0 γ * γ * (−Q 2 , −Q 2 ) ∼ 1/Q 4 . The LMD model is constructed in such a way that it fulfills the OPE constraint, but fails to reproduce the Brodsky-Lepage behavior. Finally, the LMD+V models is a generalization of the LMD model and contains two vector resonances ρ and ρ . It can be made to fulfill both the BL and the OPE constraints, for the prize of a large number of free parameters.
In order to reduce the number of fit parameters for all the models, a global fit is performed where all lattice ensembles are fitted simultaneously assuming a linear dependence of each model parameter on the lattice spacing a and on the squared pion mass m 2 π , see Ref. [21] for all the details and results of the fits.
The fits for the VMD and the LMD model are also used to perform the integration in Eq. (14) up to infinite τ, see Fig. 4 (right). For |τ| < 1.3 fm the lattice data are used. The dependence on these models for large τ is small, but the behavior for small τ is very different. In fact, both the lattice data and the LMD model show a cusp at τ = 0, which is related to the OPE in the double-virtual case, whereas the VMD model is smooth at τ = 0, since the VMD model falls off too fast at large Q 2 .
The VMD model leads to a poor description of our data, with χ 2 /d.o.f. = 2.9 (uncorrelated fit), especially in the double-virtual case and at large Euclidean momenta, see Fig. 5. The normalization F VMD π 0 γ * γ * (0, 0) is off from the value expected from the chiral anomaly and the fitted vector meson mass does not agree with the ρ-mass.
On the other hand, both the LMD model and the LMD+V model lead to a quite good fit. The LMD model has a χ 2 /d.o.f. = 1.3 (uncorrelated fit) and leads to a determination of the chiral anomaly at the 7% level, with a value in agreement with expectations. The fitted vector meson mass is close to the rho-meson mass. Furthermore, the fit result for another parameter that is related to the OPE is compatible with the theoretical expections. Despite the fact that the LMD model fails to reproduce the BL behavior for the single-virtual form factor, this does not seem to affect the global fit, since there are only few lattice data points at rather low momenta in the single virtual case, see Figs. 4 and 5, i.e. one is not yet sensitive to the asymptotic behavior.
Finally, the LMD+V model has a χ 2 /d.o.f. = 1.4 (uncorrelated fit), however, only after the model parameters that are related to the two vector meson masses, the BL behavior and the OPE constraint, have been fixed, otherwise no stable fit was obtained. The chiral anomaly is reproduced with 9% accuracy and two other fitted model parameters are close to results obtained in phenomenological analyses of the TFF [33].
The form factors in the three models, extrapolated to the physical point, are shown in Fig. 6. In the single-virtual case, the LMD+V model is in quite good agreement with the experimental data from Ref. [28]. The LMD model starts to deviate already at Q 2 = 1 GeV 2 . In the doublevirtual case, the LMD and LMD+V models are quite similar and already close to their asymptotic behavior at the largest point Q 2 ∼ 1.5 GeV 2 where we have lattice data. Using the LMD+V model at the physical point from our fit, the pion-pole contribution to HLbL in the muon g − 2 can be obtained from the three-dimensional integral representation from Ref. [1], with the result [21], Note that the given error is only statistical. No attempt has been made to estimate the systematical errors from using different fit models. Since the relevant momentum range in the pion-pole contribution is below 1 GeV  [GeV] [GeV] Figure 6. Extrapolation of the lattice data to the continuum and the physical pion mass for the VMD, LMD and LMD+V models.
(Left) Single-virtual form factor, compared with experimental results [28] and the expectation for the asymptotic behavior according to Brodsky-Lepage [30]. (Right) Double-virtual form factor at Q 2 1 = Q 2 2 and the expectation from the OPE [31].

HLbL forward scattering amplitudes in lattice QCD
Using parity and time-reversal invariance of QCD, there are eight independent light-by-light forward scattering amplitudes describing the process γ * (λ 1 , q 1 ) γ * (λ 2 , q 2 ) → γ * (λ 1 , q 1 ) γ * (λ 2 , q 2 ), where q i and λ ( ) i = 0, ± are the momenta and helicities of the virtual photons (i = 1, 2). Six amplitudes are even and two are odd functions of the crossing-symmetric variable ν = q 1 · q 2 . The forward scattering amplitudes M λ 1 λ 2 λ 1 λ 2 can be related via the optical theorem to two-photon fusion amplitudes M λ 1 λ 2 for the process γ * (λ 1 , q 1 ) γ * (λ 2 , q 2 ) → X(p X ) as follows: Unitarity and analyticity then allow one to write dispersion relations in ν at fixed values of the virtualities Q 2 i = −q 2 i . Performing one subtraction to get a faster convergence and to suppress higher resonance states, one obtains the following sum rules (omitting all dependence on the virtualities and the helicities): For later use, we introduce the following notation for the subtracted even M(q 2 1 , q 2 2 , ν) = M(q 2 1 , q 2 2 , ν) − M(q 2 1 , q 2 2 , 0) and odd M(q 2 1 , q 2 2 , ν) = M(q 2 1 , q 2 2 , ν) − νM (q 2 1 , q 2 2 , 0) amplitudes. As proposed first in Ref. [23] for the forward scattering amplitude M T T and extended recently to all eight amplitudes in Ref. [24], the scattering amplitudes on the lefthand sides of Eqs. (17)-(18) can be computed from the correlation function of four vector currents on the lattice. On the other hand, the right-hand sides of Eqs. (17)- (18) are related to the two-photon fusion processes using Eq. (16). The latter are described by single-meson TFFs which can be parametrized by simple models, as sketched below. Fitting the lattice data then determines the model parameters and the TFFs. This will allow one to calculate the contributions from single-meson poles to the HLbL contribution in the muon g − 2.
The lattice simulations are performed on five CLS lattice ensembles [27] with two degenerate light dynamical quarks at two lattice spacings a = 0.048 fm (one ensemble) and a = 0.065 fm (four ensembles) and pion masses in the range from 194−437 MeV. For all ensembles the fully connected and for two ensembles also the leading quarkdisconnected diagrams contributing to the four-point function are taken into account. For each ensemble, the correlation function is computed at up to three values of Q 2 1 below 0.8 GeV 2 and for all values of Q 2 2 ≤ 4 GeV 2 . The result for a subset of four scattering amplitudes is shown in Fig. 7.
The photon fusion reaction on the right-hand sides of Eqs. (17)- (18) can produce any C-parity-even state X. The main contribution is expected from the pseudoscalar (0 −+ ), scalar (0 ++ ), axial-vector (1 ++ ) and tensor (2 ++ ) mesons, where we consider in each channel only the lightest state. We are working with two degenerate dynamical quarks and fit our phenomenological model to only the fullyconnected diagrams. To compensate for this, we include only the contributions from isovector mesons, multiplied by a factor of 34/9 [10,24].
There is one TFF for the pseudoscalars, two each for the scalars and axial-vectors and four for the tensor mesons [22]. For the pseudoscalars, we take the TFF from Ref. [21] evaluated on the same lattice ensembles. The other TFFs are parametrized as follows where we assume a monopole ansatz (n = 1) for the scalars and a dipole ansatz (n = 2) for the axial-vectors and tensor mesons, parametrized by the mass scale Λ X . We assume one common mass for the scalar TFFs, one mass for the TFFs of the axial-vectors and four different masses for the TFFs of the tensor mesons. These six masses will be considered as free fit parameters. The normalization of the TFFs is given by the twophoton decay width and is taken from experiment where available, e.g. for the scalars Γ γγ = πα 2 4 m S F T Sγ * γ * (0, 0) 2 (an appropriately defined effective two-photon width is employed for the axial-vectors [22]). Since not all normalizations have been measured, further input from dispersive sum rules from Ref. [15] is used. Furthermore, for the pseudoscalars, again the lattice data from Ref. [21] are used.
As an example, the pseudoscalar contribution to the cross-section σ T T of two transversely polarized photons is given in the narrow-width approximation by where X = ν 2 − Q 2 1 Q 2 2 is the virtual-photon flux factor. Similar results can be obtained for the other mesons where we assume a Breit-Wigner shape for the resonances [24].
In Fig. 7 the results for four scattering amplitudes of a combined fit to all eight amplitudes with the phenomenological model described above is shown for one lattice ensemble. For four ensembles, the χ 2 /d.o.f. is quite good, between 1.13 − 1.35. The fit for the ensemble with the heaviest pion mass m π = 437 MeV is, however, not very satisfactory and that ensemble is left out in the chiral extrapolation to the physical pion mass.
The relative contribution of the different mesons to the individual scattering amplitudes was also studied in Ref. [24]. The pseudoscalar and tensor mesons give the dominant contribution to the amplitudes M T T , M   the main contribution are from axial and tensor mesons. In the amplitudes M a T L , M τ T L , M LL scalar, axial and tensor mesons contribute significantly. The contribution from γ * γ * → π + π − , evaluated with scalar QED dressed with a monopole vector form factor, is always small compared to the other channels.
The lattice simulations are performed for ensembles away from the physical quark masses. The pion and ρmeson masses are set to their lattice values obtained from the exponential decay of the pseudoscalar and vector twopoint functions. For other resonances, we assume a constant shift in the masses m X = m phys X + (m lat ρ − m exp ρ ). In Table 2 we compare the results of the masses in the TFFs, extrapolated to the physical pion mass, to experimental and phenomenological determinations [15,35]. The agreement is reasonably good for M S , M (2) T , M (0,L) T , although the scalar mass on the lattice is a bit high. There are quite strong tensions for M A , M (1) T , M (0,T ) T . In particular the latter two tensor masses on the lattice are almost a factor two larger than the phenomenological determinations. See Ref. [24] for a more detailed discussion and potential reasons for the disagreements. Note in particular, that we have not yet performed the continuum limit.
Overall, we get a good description of the lattice data with the lattice determination of the pion TFF [21] and the simple monopole or dipole ansätze from Eq. (19) for the various TFF's with one resonance in each channel. Table 2. Chiral extrapolation to the physical pion mass for the scalar monopole mass M S , the axial dipole mass M A and the four tensor dipole masses corresponding to different helicities, compared to experimental or phenomenological determinations. All masses are given in GeV. The ensembles have a finite lattice spacing of a = 0.065 fm and the ensemble with the largest pion mass has been excluded in the chiral extrapolation.   (19) 0.877(66)

Conclusions
Lattice QCD can provide a model-independent, firstprinciple calculation of the HLbL contribution to the muon g − 2. Recent first preliminary results by RBC-UKQCD [17] and the lattice group at Mainz [18-21, 23, 24] look very promising. Hopefully, in a few years time when the final result from the Fermilab experiment will be published, an estimate with 10% uncertainty (combined statistical and controlled systematics errors) can be reached, which would match the expected experimental precision, if one assumes that the central value for HLbL is close to current model estimates a HLbL µ ≈ 100 × 10 −11 [1,2,11,12].
Since HLbL involves a rank-four tensor with many independent momenta (or space-time points in position space), it is a very complicated object. We therefore think it makes sense to have as many tests and observables as possible on the lattice, not just the final number for a HLbL µ itself. Therefore the calculation of the pion (light pseudoscalar) transition form factor and of HLbL forward scattering amplitudes, combined with HLbL sum rules, will be a valuable tool to compare different lattice calculations, once they become available.