Hadronic light-by-light scattering contribution to the muon 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 , 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 in the spacelike region with photon virtualities up to 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 fusion cross sections and then compared with phenomenological models.
The anomalous magnetic moment of the muon has served for many years as a precision test of the Standard Model [1, 2, 4, 3] and it has also played an important role in many presentations at this meeting. There is a discrepancy of 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 , 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 . But also lattice QCD is getting more and more precise , and, hopefully, also a new method using muon-electron scattering to measure the running of and the HVP in the spacelike region , will be feasible at some point with the required precision.
On the other hand, the HLbL contribution to the muon , see Fig. 1, has only been calculated using models so far [1, 4, 10] and the frequently used estimates from Refs. [11, 12, 1] (revised slightly in Ref. ) 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. ) to determine the presumably numerically dominant contributions from a single neutral pion-pole (light pseudoscalar-pole) and from the two-pion intermediate state (pion-loop) based on input from experimental data for the dispersion 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 and some promising progress has been achieved recently by the RBC-UKQCD collaboration .
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 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  and we analyzed HLbL forward scattering amplitudes that can be compared using dispersion relations (HLbL sum rules)  to phenomenological models for photon fusion processes [23, 18, 24]. These three topics will be discussed in the following three sections. For all the details, we refer to the quoted papers.
2 Position-space approach to HLbL in the muon 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 kernel that describes the muon and photon propagators and a non-perturbative QCD four-point function (denoted by a blob in the Feynman diagram) that will be evaluated on the lattice.
with the spatial moment of the four-point function
We evaluate the QED kernel in the continuum and in infinite volume and thereby avoid 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 and .
The QED kernel can be decomposed into several tensors
The are traces of gamma matrices and just yield sums of products of Kronecker deltas. The tensors are decomposed into a scalar , vector and tensor part
These can be parametrized by six weight functions
that depend on the three variables , and . 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 has been given in Ref. . In Fig. 2 we show the two weight functions and as a function of , for a fixed value of and three values of . Plots for all six weight functions can be found in Ref. .
To test our semi-analytical expressions and the software for the QED kernel, we have computed the -pole contribution to HLbL in a simple vector-meson dominance (VMD) model as well as the lepton-loop contribution to in QED, where the results are well known. To this aim, we first derived analytical expressions for these contributions to the four-point function in position-space, see Ref.  for details.
In figure 3 we plot the integrand of the final integration over of the HLbL contribution (LbL in QED) for these two examples, after contracting all Lorentz indices and the integrations over and have been performed in Eq. (1).
From the plot of the integrand for the pion-pole contribution in Fig. 3 one observes that this contribution to is remarkably long-range with a long negative tail at large . One expects an exponential decay of the correlation function. But this seems to be countered by some non-negligible power-like behavior . For pion masses we reproduce the known results, which can be easily obtained from the three-dimensional integral representation in momentum space given in Ref. , at the percent level. On the other hand, for the physical pion mass, one will need rather large lattices of the order of to capture the negative tail at large 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.  and Section 3.
For the lepton-loop in QED the behavior of the integrand for small is compatible with . 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 of the lepton in the loop. Furthermore there is again a long negative tail at large 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 we reproduce the analytically known results for in QED  at the percent level, see Table 1. On the other hand, for the lightest lepton mass , some further refinements of our numerical evaluation are needed. Once this is achieved, we plan to make the QED kernel publicly available.
3 Lattice calculation of the pion transition form factor
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. 
provided the photon virtualities satisfy 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 denotes the zeroth (energy) component of the four-momentum .
The main object to compute on the lattice is the three-point function
Here is the time separation between the two vector currents and . The matrix element in Eq. (10) with an on-shell pion is obtained by considering the limit of large . With the definitions
The overlap factor and the pion energy can be obtained from the asymptotic behavior of the pseudoscalar two-point function.
For the calculation of on the lattice in Ref.  we used eight CLS (Coordinated Lattice Simulations) lattice ensembles  with dynamical quarks with three different lattice spacings and pion masses in the range .
We choose the pion rest frame , with photons back-to-back spatially , where the kinematical range accessible on the lattice is given by , . With multiple values of , we obtain mostly spacelike photon virtualities up to , as can be seen in Fig. 4 (left). In practice, discrete values of have been used to sample the momenta. Note that on the lattice it is actually easier to access the double-virtual TFF , in particular near the diagonal , than the single-virtual form factor along the two axis, in contrast to the situation in experiments .
From the theoretical side, the form factor is constrained by the chiral anomaly such that  (in the chiral limit). For the single-virtual form factor one expects the Brodsky-Lepage (BL) behavior . 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 in the continuum and for the physical pion mass, we fit our lattice data, obtained for the eight ensembles with different lattice spacings and pion masses , 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 single-virtual case (and describes quite well the available experimental data below ), but falls off too fast for the double virtual case . 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 and on the squared pion mass , see Ref.  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 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 , which is related to the OPE in the double-virtual case, whereas the VMD model is smooth at , since the VMD model falls off too fast at large .
The VMD model leads to a poor description of our data, with (uncorrelated fit), especially in the double-virtual case and at large Euclidean momenta, see Fig. 5. The normalization 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 (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 (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 .
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. . The LMD model starts to deviate already at . In the double-virtual case, the LMD and LMD+V models are quite similar and already close to their asymptotic behavior at the largest point 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 can be obtained from the three-dimensional integral representation from Ref. , with the result ,
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 , i.e. where most of our lattice data points are obtained, the extrapolation to large momenta with the LMD+V or the LMD model does not have a large effect on the final result. The result in Eq. (15) is fully consistent with most model calculations which yield results in the range , but with rather arbitrary, model-dependent error estimates, see Refs. [1, 10, 34] and references therein.
4 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 , where and are the momenta and helicities of the virtual photons . Six amplitudes are even and two are odd functions of the crossing-symmetric variable . The forward scattering amplitudes can be related via the optical theorem to two-photon fusion amplitudes for the process as follows:
Unitarity and analyticity then allow one to write dispersion relations in at fixed values of the virtualities . 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):
with . For later use, we introduce the following notation for the subtracted even and odd amplitudes.
As proposed first in Ref.  for the forward scattering amplitude and extended recently to all eight amplitudes in Ref. , the scattering amplitudes on the left-hand 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 .
The lattice simulations are performed on five CLS lattice ensembles  with two degenerate light dynamical quarks at two lattice spacings (one ensemble) and (four ensembles) and pion masses in the range from . For all ensembles the fully connected and for two ensembles also the leading quark-disconnected 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 below and for all values of . 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 . The main contribution is expected from the pseudoscalar (), scalar (), axial-vector () and tensor ( 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 fully-connected diagrams. To compensate for this, we include only the contributions from isovector mesons, multiplied by a factor of [10, 24].
There is one TFF for the pseudoscalars, two each for the scalars and axial-vectors and four for the tensor mesons . For the pseudoscalars, we take the TFF from Ref.  evaluated on the same lattice ensembles. The other TFFs are parametrized as follows
where we assume a monopole ansatz for the scalars and a dipole ansatz for the axial-vectors and tensor mesons, parametrized by the mass scale . 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 two-photon decay width and is taken from experiment where available, e.g. for the scalars (an appropriately defined effective two-photon width is employed for the axial-vectors ). Since not all normalizations have been measured, further input from dispersive sum rules from Ref.  is used. Furthermore, for the pseudoscalars, again the lattice data from Ref.  are used.
As an example, the pseudoscalar contribution to the cross-section of two transversely polarized photons is given in the narrow-width approximation by
where 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 .
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 is quite good, between . The fit for the ensemble with the heaviest pion mass 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. . The pseudoscalar and tensor mesons give the dominant contribution to the amplitudes , and that involve two transverse photons. The pseudoscalar meson does not contribute to , where the main contribution are from axial and tensor mesons. In the amplitudes , , 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 two-point functions. For other resonances, we assume a constant shift in the masses . In Table 2 we compare the results of the masses in the TFFs, extrapolated to the physical pion mass, to experimental and phenomenological determinations [35, 15]. The agreement is reasonably good for , , although the scalar mass on the lattice is a bit high. There are quite strong tensions for . In particular the latter two tensor masses on the lattice are almost a factor two larger than the phenomenological determinations. See Ref.  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  and the simple monopole or dipole ansätze from Eq. (19) for the various TFF’s with one resonance in each channel.
Lattice QCD can provide a model-independent, first-principle calculation of the HLbL contribution to the muon . Recent first preliminary results by RBC-UKQCD  and the lattice group at Mainz [18, 19, 20, 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 [11, 12, 1, 2].
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 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.
We are grateful for the use of CLS lattice ensembles. We acknowledge the use of computing time on the JUGENE and JUQUEEN machines located at Forschungszentrum Jülich, Germany. The correlation functions were computed on the clusters “Wilson” at the Institute of Nuclear Physics, University of Mainz, “Clover” at the Helmholtz-Institute Mainz and “Mogon” at the University of Mainz, Germany. This work was partially supported by the Deutsche Forschungsgemeinschaft (DFG) through the Collaborative Research Center “The Low-Energy Frontier of the Standard Model” (SFB 1044).
-  F. Jegerlehner and A. Nyffeler, Phys. Rept. 477, 1 (2009).
-  F. Jegerlehner, EPJ Web Conf. 118, 01016 (2016); F. Jegerlehner, arXiv:1705.00263 [hep-ph].
-  S. Laporta, Phys. Lett. B 772, 232 (2017); S. Laporta, talk at this meeting.
-  M. Knecht, talk at this meeting.
-  D. Stöckinger, talk at this meeting.
-  T. Gorringe, G. Marshall, talks at this meeting.
-  D. Nomura, V. Sauli, F. Ignatov, A. Denig, talks at this meeting.
-  M. Marinkovic, C. Lehner, talks at this meeting.
-  L. Trentadue, U. Marconi, P. Mastrolia, F. Piccinini, talks at this meeting.
-  J. Bijnens and J. Relefors, JHEP 1609, 113 (2016); J. Bijnens, talk at this meeting, arXiv:1712.09787 [hep-ph].
-  J. Prades, E. de Rafael and A. Vainshtein, Adv. Ser. Direct. High Energy Phys. 20, 303 (2009).
-  A. Nyffeler, Phys. Rev. D 79, 073012 (2009).
-  G. Colangelo et al., JHEP 1409, 091 (2014); G. Colangelo et al., Phys. Lett. B 738, 6 (2014); G. Colangelo et al., JHEP 1509, 074 (2015); Phys. Rev. Lett. 118, 232001 (2017); JHEP 1704, 161 (2017); G. Colangelo, talk at this meeting.
-  V. Pauk and M. Vanderhaeghen, arXiv:1403.7503 [hep-ph]; Phys. Rev. D 90, 113012 (2014).
-  I. Danilkin and M. Vanderhaeghen, Phys. Rev. D 95, 014019 (2017).
-  F. Hagelstein and V. Pascalutsa, arXiv:1710.04571 [hep-ph].
-  M. Hayakawa, et al., PoS LAT 2005, 353 (2006); T. Blum, M. Hayakawa and T. Izubuchi, PoS LATTICE 2012, 022 (2012); T. Blum et al., Phys. Rev. Lett. 114, 012001 (2015); T. Blum et al., Phys. Rev. D 93, 014503 (2016); T. Blum et al., Phys. Rev. Lett. 118, 022005 (2017); Phys. Rev. D 96, 034515 (2017); C. Lehner, talk at this meeting.
-  J. Green et al., PoS LATTICE 2015, 109 (2016).
-  N. Asmussen et al., PoS LATTICE 2016, 164 (2016).
-  N. Asmussen et al., arXiv:1711.02466 [hep-lat].
-  A. Gérardin, H. B. Meyer and A. Nyffeler, Phys. Rev. D 94, 074507 (2016).
-  V. Pascalutsa and M. Vanderhaeghen, Phys. Rev. Lett. 105, 201603 (2010); V. Pascalutsa, V. Pauk and M. Vanderhaeghen, Phys. Rev. D 85, 116001 (2012).
-  J. Green et al., Phys. Rev. Lett. 115, 222003 (2015).
-  A. Gérardin et al., arXiv:1712.00421 [hep-lat].
-  S. Laporta and E. Remiddi, Phys. Lett. B 301, 440 (1993); M. Passera, private communication.
-  X. d. Ji and C. w. Jung, Phys. Rev. Lett. 86, 208 (2001); Phys. Rev. D 64 (2001) 034506; J. J. Dudek and R. G. Edwards, Phys. Rev. Lett. 97, 172001 (2006); X. Feng et al., Phys. Rev. Lett. 109, 182001 (2012).
-  P. Fritzsch et al., Nucl. Phys. B 865, 397 (2012).
-  H. J. Behrend et al. [CELLO Collaboration], Z. Phys. C 49, 401 (1991); J. Gronberg et al. [CLEO Collaboration], Phys. Rev. D 57, 33 (1998); B. Aubert et al. [BaBar Collaboration], Phys. Rev. D 80, 052002 (2009); S. Uehara et al. [Belle Collaboration], Phys. Rev. D 86, 092007 (2012).
-  S. L. Adler, Phys. Rev. 177, 2426 (1969); J. S. Bell and R. Jackiw, Nuovo Cim. A 60, 47 (1969).
-  G. P. Lepage and S. J. Brodsky, Phys. Lett. B 87, 359 (1979); Phys. Rev. D 22, 2157 (1980); S. J. Brodsky and G. P. Lepage, Phys. Rev. D 24, 1808 (1981).
-  V. A. Nesterenko and A. V. Radyushkin, Sov. J. Nucl. Phys. 38, 284 (1983); V. A. Novikov et al., Nucl. Phys. B 237, 525 (1984).
-  B. Moussallam, Phys. Rev. D 51, 4939 (1995); M. Knecht et al., Phys. Rev. Lett. 83, 5230 (1999).
-  M. Knecht and A. Nyffeler, Eur. Phys. J. C 21, 659 (2001); K. Melnikov and A. Vainshtein, Phys. Rev. D 70, 113006 (2004).
-  M. Knecht and A. Nyffeler, Phys. Rev. D 65, 073034 (2002); A. Nyffeler, Phys. Rev. D 94, 053006 (2016).
-  P. Achard et al. [L3 Collaboration], Phys. Lett. B 526, 269 (2002); P. Achard et al. [L3 Collaboration], JHEP 0703, 018 (2007); M. Masuda et al. [Belle Collaboration], Phys. Rev. D 93, 032003 (2016).