Hadronic light-by-light scattering amplitudes from lattice QCD versus dispersive sum rules

The hadronic contribution to the eight forward amplitudes of light-by-light scattering ($\gamma^*\gamma^*\to \gamma^*\gamma^*$) is computed in lattice QCD. Via dispersive sum rules, the amplitudes are compared to a model of the $\gamma^*\gamma^*\to {\rm hadrons}$ cross sections in which the fusion process is described by hadronic resonances. Our results thus provide an important test for the model estimates of hadronic light-by-light scattering in the anomalous magnetic moment of the muon, $a_\mu^{\rm HLbL}$. Using simple parametrizations of the resonance $M\to \gamma^*\gamma^*$ transition form factors, we determine the corresponding monopole and dipole masses by performing a global fit to all eight amplitudes. Together with a previous dedicated calculation of the $\pi^0\to \gamma^*\gamma^*$ transition form factor, our calculation provides valuable information for phenomenological estimates of $a_\mu^{\rm HLbL}$. The presented calculations are performed in two-flavor QCD with pion masses extending down to 190\,MeV at two different lattice spacings. In addition to the fully connected Wick contractions, on two lattice ensembles we also compute the (2+2) disconnected class of diagrams, and find that their overall size is compatible with a parameter-free, large-$N$ inspired prediction, where $N$ is the number of colors. Motivated by this observation, we estimate in the same way the disconnected contribution to $a_\mu^{\rm HLbL}$.

The hadronic contribution to the eight forward amplitudes of light-by-light scattering (γ * γ * → γ * γ * ) is computed in lattice QCD. Via dispersive sum rules, the amplitudes are compared to a model of the γ * γ * → hadrons cross sections in which the fusion process is described by hadronic resonances. Our results thus provide an important test for the model estimates of hadronic light-by-light scattering in the anomalous magnetic moment of the muon, a HLbL µ . Using simple parametrizations of the resonance M → γ * γ * transition form factors, we determine the corresponding monopole and dipole masses by performing a global fit to all eight amplitudes. Together with a previous dedicated calculation of the π 0 → γ * γ * transition form factor, our calculation provides valuable information for phenomenological estimates of a HLbL µ . The presented calculations are performed in two-flavor QCD with pion masses extending down to 190 MeV at two different lattice spacings. In addition to the fully connected Wick contractions, on two lattice ensembles we also compute the (2+2) disconnected class of diagrams, and find that their overall size is compatible with a parameter-free, large-N inspired prediction, where N is the number of colors. Motivated by this observation, we estimate in the same way the disconnected contribution to a HLbL The non-vanishing probability of two photons scattering off each other is a striking prediction of quantum electrodynamics (QED) [1,2]. The smallness of the cross section has so far prohibited a direct experimental observation, although evidence for the phenomenon has recently been found by the ATLAS experiment in relativistic heavy-ion collisions [3]. Equally interesting is the scattering of spacelike virtual photons, γ * γ * → γ * γ * . While the contributions of virtual leptons are calculable in QED perturbatively, the hadronic contributions require a nonperturbative approach. When the photons are real or spacelike, dispersive sum rules can be used to express the forward γ * γ * scattering amplitudes in terms of experimentally more accessible γ * γ * -fusion cross sections [4][5][6][7]. The hadronic contributions to the γ * γ * → γ * γ * amplitudes can also be computed ab initio using lattice QCD [8].
One very timely application of hadronic light-by-light (HLbL) scattering is the anomalous magnetic moment of the muon, a µ = 1 2 (g − 2) µ . The current discrepancy between the direct measurement of a µ and the Standard Model prediction amounts to about 3.6 standard deviations [9]. While the current theory and experimental errors are comparable in size, two new (g − 2) µ experiments [10,11] in preparation at Fermilab and J-PARC aim at reducing the experimental error by a factor of four. The largest sources of theory error are contributions from the hadronic vacuum polarization (HVP) and from HLbL scattering. The latter is expected to dominate in the future in view of the dedicated measurements at e + e − colliders ever better constraining the former. Although an active area of research, the experimental data needed for the recently proposed data-driven dispersive approaches to HLBL [12][13][14][15][16] are harder to obtain, and lattice QCD calculations are in particularly high demand.
Lattice QCD calculations of hadron structure have been steadily advancing in recent years. Several collaborations calculate hadronic observables directly at physical values of the quark masses. At least two collaborations are addressing the HLbL contribution to (g − 2) µ on the lattice [8,[17][18][19][20][21][22]. Although the calculation poses serious challenges due to the complexity of the four-point function and the long-range nature of the dominant contribution, as a spacelike quantity it is well suited for a first-principles treatment directly in the Euclidean theory. A second role of lattice QCD is to provide the necessary hadronic input for the model and dispersive approaches to a HLbL µ . Model calculations (see e.g. [23] for a recent overview) consistently suggest that the contribution of the pseudoscalar mesons (π 0 , η, η ) is dominant, and therefore determining their respective transition form factors is of primary importance. A first calculation of the π 0 → γ * γ * transition form factor in the range of photon virtualities relevant to (g − 2) µ has been carried out on the lattice [24]. An extension of this calculation to the η, η mesons is possible, though more demanding, due to the appearance of disconnected Wick-contraction diagrams. Computing the spectrum and two-photon coupling of the scalar, axial-vector and tensor mesons is qualitatively more complicated in lattice QCD, since these states are resonances and require a dedicated treatment.
In many ways, the light-by-light scattering amplitudes are the more accessible observable in lattice calculations, because they involve spacelike photons that can be treated directly in the Euclidean theory. The lattice calculation of the cross section γ * γ * → ππ is for instance more complex than calculating the cross section γ * γ * → γ * γ * for spacelike photons. In experiments, the weakness of the electromagnetic coupling would make such a measurement impractical, but in lattice QCD the factor e 4 merely multiplies a four-point correlation function at the end of the calculation.
In this article, we compute the HLbL scattering amplitude for spacelike photons in lattice QCD. Being parametrized by functions of six Lorentz invariants, it is a complicated object. We focus on the forward amplitudes because they are simpler functions of three invariants and, using the optical theorem, they are related to γ * γ * → hadrons cross sections. Our objectives are: (i) Provide a stringent test that the light-by-light amplitude for spacelike photons is correctly described by the type of hadronic model used so far to estimate a HLbL µ . The model includes the exchange of pseudoscalar, scalar, axial-vector and tensor mesons.
(ii) Provide information on their two-photon transition form factors via global fits.
(iii) Compare the transition form factors to phenomenological determinations based on lightby-light sum rules.
In [8], we laid out the method and computed the forward amplitude sensitive to the total transverse γ * γ * cross section, M TT , via a dispersive sum rule. Here we extend the comparison between lattice data and phenomenological parametrizations of the γ * γ * → hadrons cross sections to encompass all eight forward amplitudes. This more extensive analysis allows us to place much stronger constraints on the size of the contributions of different resonances, because they contribute to different amplitudes with different weight factors, often even with opposite signs. Our lattice calculation is performed in QCD with two flavors of light quarks; it involves pion masses down to 190 MeV and two lattice spacings. While a fully realistic lattice calculation would have to include at least a dynamical strange quark, the present calculation does provide a suitable test of hadronic models via dispersive sum rules, since at the required level of precision it is fairly straightforward to adapt these models to QCD without the strange quark, as discussed in section V.
Here we have computed the fully connected and the dominant disconnected Wick-contraction diagrams. We discuss what these classes of diagrams correspond to in terms of the quantum numbers of the exchanged resonances. The large-N inspired approximation that a quark loop containing a single, vector-current insertion gives a negligible contribution, corresponds, in twoflavor QCD, to including only isovector resonances, enhanced by a factor of 34/9. This interpretation of the fully connected class of diagrams was first pointed out in [25,26], mainly concerning the pseudoscalar sector; see also the arguments presented in [27]. We rederive the result in detail in the SU(2) and in the SU(3) flavor symmetric theory under slightly weaker assumptions; see section III B.
The paper is organized as follows. We begin by introducing the theoretical background for light-by-light scattering in section II. We then describe the lattice method for computing the hadronic light-by-light amplitude, including the analytic continuation and the numerical method to obtain the fully connected and the (2+2) disconnected four-point function (section III). Section IV presents our numerical results in two-flavor QCD, while some additional material is provided in appendix B. After introducing the details of the hadronic model for the γ * γ * → hadrons cross sections in section V and appendix A, we perform fits to the lattice data in section VI. We compare the results for the transition form factors to existing phenomenological estimates. In section VII, we discuss our results on the leading disconnected diagrams to the HLbL amplitude and present what results would have to be obtained in lattice QCD for the connected and the leading disconnected diagram contributions to a HLbL µ in order to confirm the state-of-the-art model estimate. We conclude in section VIII.

II. FORWARD LIGHT-BY-LIGHT SCATTERING AND SUM RULES
In order to establish our notation, we start by recalling the dispersive sum rules for the scattering of spacelike photons [4,5]. Just as for real photons [28], they are based on unitarity and analyticity of the forward scattering amplitude. More specifically, the optical theorem allows one to relate the absorptive part of the γ * (λ 1 , q 1 )γ * (λ 2 , q 2 ) → γ * (λ 1 , q 1 )γ * (λ 2 , q 2 ) forward scattering amplitude to fusion cross sections for the process γ * γ * → X, where X stands for any C-parity even final state. The relevant kinematic variables are the photon virtualities, q 2 i = −Q 2 i , (i = 1, 2), and the crossing-symmetric variable ν = q 1 · q 2 , which is related to the squared centerof-mass energy by s = 2ν − Q 2 1 − Q 2 2 . Denoting the absorptive part of the helicity amplitude M λ 1 λ 2 ,λ 1 λ 2 by The largest value of |ν| that can be reached with Euclidean kinematics is limited by the virtualities of the photons 3 , |ν| ≤ (Q 2 1 Q 2 2 ) 1/2 ≤ 1 2 (Q 2 1 + Q 2 2 ) ≡ ν 0 , while the nearest singularity is the s-channel π 0 pole located at ν π = 1 2 (m 2 π + Q 2 1 + Q 2 2 ). 1. Special case of (anti)parallel momenta The tensors T E , resulting from Eq. (5) translated to Euclidean space using Eq. (15), are not defined for collinear Q 1 and Q 2 , since in that case X = 0. However, if we start with non-collinear momenta and rotate Q 2 toward being (anti)parallel with Q 1 , then each tensor has such a limit. This limit depends on the initial direction of Q 2 ; we will use the average over this direction to define T E in the collinear limit.
We define the projector Here V is a unit vector orthogonal to Q 1 , pointing in the direction from which Q 2 approached being collinear with Q 1 . Since K 2 · Q 1 = i −X/Q 2 2 → 0, in the collinear limit any contraction of T E µνµ ν with Q 1 will vanish. Thus, after averaging over all V orthogonal to Q 1 , each T E µνµ ν will be a linear combination of We obtain the prefactors by contracting the indices in three different ways. For this, we will make use of Denoting by . . . V the average over V , we find B. Flavor structure of the four-point function In numerical lattice QCD calculations of n-point functions, the quark path integral is evaluated analytically to yield a sum of contractions of quark propagators. For the four-point function of vector currents, these fall into five distinct topologies, illustrated in Fig. 1.
The calculation of all Wick-contraction topologies is demanding. In many instances, disconnected diagram contributions have been found to make numerically small contributions to hadronic matrix elements, though not always [32]. Quark loops generated by a single vector current have been empirically found to be particularly suppressed (see for instance [24,[33][34][35]). At short distances, perturbation theory provides an explanation for the suppression of this type of contribution, since it requires the exchange of at least three gluons [36]. On the other hand, it is well known that the disconnected diagram is responsible for the difference between the pion and the η mass in the pseudoscalar two-point function, and is therefore crucial at long distances.
The importance of the disconnected diagrams in the HLbL amplitude has been pointed out in [25,26], showing that the pion and η pole contributions would have the wrong weight factors if only the connected diagrams were included. Here, 1. we use (a) flavor symmetry, either SU(2) or SU(3), and (b) the assumption that Wickcontraction diagrams where a vector current appears as the only insertion in a quark loop, thus producing a factor Tr{γ µ S(X, X)}, are negligible; 2. we then derive the weight factors of non-singlet and singlet mesons in the fully connected and the (2+2) disconnected contribution to the HLbL amplitude; 3. we show that whenever the HLbL amplitude is dominated by the pole-exchange of an isovector resonance, isospin symmetry induces relations between different Wick-diagram topologies.
The main result is that under the assumptions stated under (1.), in the fully connected diagrams the non-singlet meson poles over-contribute by a factor 34/9 (respectively a factor 3) in QCD with two (respectively three) degenerate flavors of quarks, while the singlet mesons do not contribute. The (2+2) disconnected diagrams contain the singlet-meson contribution and correct the fully connected diagram by compensating with (−25/9) (respectively −2) times the non-singlet meson contribution. For QCD with a realistic quark spectrum, we expect the relations between the classes of diagrams to lie between the quoted predictions.
The starting point is the observation, based on Fig. 1, that the Wick contractions contributing to the HLbL amplitude can be written as (we drop the space-time arguments and indices of the four-point amplitudes) In the case that the quark masses are all equal, one can drop the flavor indices, e.g. Π 2+2 f 1 ,f 2 → Π 2+2 . In particular, let Π be the four-point function of the 'up' currentūγ µ u. Since in that case f Q n f = 1 ∀n, and because all quark lines carry the same quark mass, the Π X f are precisely the Wick contractions appearing in Π (X = 4, 2 + 2, 3 + 1, 2 + 1 + 1, 1 + 1 + 1 + 1), including the normalization. On the other hand, the HLbL amplitude is expressed as Wick contractions Π X f weighted by polynomials in the quark charges. An important observation is that with two flavors of quarks, there are only three linearly independent symmetric polynomials in Q u and Q d , for instance With three flavors, there are four such independent polynomials in (Q u , Q d , Q s ), for instance We assume exact isospin SU(2) symmetry. The photon couples to the electromagnetic current, whose isospin decomposition reads J e.m.
The upper index on the current indicates the isospin quantum number I = 0, 1. There are both isoscalar and isovector resonances that couple to two photons. As is well known, the coupling of an isovector resonance occurs only when one of the photons couples via the isoscalar part, and one couples via the isovector part of the e.m. current 4 . Since the amplitude for a neutral pion to couple to two photons vanishes if either both are isoscalar (Q u = Q d ) or both are isovector (Q u = −Q d ), and the coupling must be quadratic in the charges, it must be proportional to (Q 2 u − Q 2 d ). Then the contributions to Π and Π HLbL of an isovector resonance M 1 are related by Correspondingly, the dependence of the transition form factor of an isoscalar resonance on the quark charges is such that it contains two independent terms, The notation indicates that F D contains all the diagrams where at least one vector current appears isolated in a quark loop. In view of the form (25) of the M 0 γγ vertex, the pole contribution of an isoscalar meson has the dependence on the quark charges, where the contributions Π B and Π C are only non-vanishing if the disconnected diagrams involving at least one isolated vector current inserted in a quark loop are non-vanishing. As discussed below, Π B is O(1/N ) and Π C is O(1/N 2 ) in the large-N power counting. By identifying the polynomials in Q u and Q d in Eqs. (24) and (20), we obtain three conditions that relate the contributions of an isovector resonance to the various Wick contractions: Thus if for a specific kinematic regime one isovector resonance exchange dominates the HLbL amplitude, it is sufficient to compute three of the five Wick-contraction classes. For instance, Π 1+1+1+1,M 1 and Π 2+1+1,M 1 can be expressed in terms of the classes Π 4,M 1 , Π 2+2,M 1 and Π 3+1,M 1 of diagrams. We also note the exact expression for the fully connected class of diagrams. Similarly, equating the expressions (26) and (20) yields the relations Eliminating Π 2+2,M 0 and Π 3+1,M 0 , we get

Large-N expectations
In terms of large-N counting, where N is the number of colors, every additional disconnected quark loop costs a factor 1/N . However it is worth looking more closely how this property emerges. For instance, just keeping the leading class of diagrams Π 4,M would not reproduce the correct dependence on Q u and Q d for the exchange of a single meson, be it isoscalar or isovector. Put in a different way, the relations (30) and (34) derived from isospin symmetry relate contributions to diagrams that scale differently with N .
The resolution of this apparent contradiction is that, at large N , neutral mesons are expected to come in degenerate pairs, one isoscalar, one isovector, due to the vanishing of the quark annihilation diagrams. Only when the sum of the contributions of a pair is considered, the large-N counting should apply. Here we only apply the large-N counting rule to the insertion of a vector current. Considering the sum of the contributions of such a pair, we obtain from (30-34) the relation Large-N counting should apply at this point, and we expect to be able to neglect in first approximation Π B , Π C (down by 1/N and 1/N 2 respectively) and the terms Π 2+1+1,M 0 +M 1 and Π 1+1+1+1,M 0 +M 1 , since they are expected to be down by 1/N 2 and 1/N 3 respectively relative to Π 4,M 0 +M 1 . All neglected terms contain at least one isolated vector current in a quark loop.
In that approximation, the leading diagram classes are given by Thus, using Eqs (36) and (24) we obtain for the fully connected contribution in Eq. (20) where we have included the physical charges of the u and d quarks in the last step. The (2+2) disconnected diagrams complement the connected diagrams to yield the full contribution, The charge factors in Eqs. (38) and (39) agree with Refs. [25,26]. The stronger large-N prediction, which however turns out not to be a good approximation in QCD, is that isoscalar and isovector states in each symmetry channel (pseudoscalars, scalars, tensors, etc.) compensate each other in the (2+2) disconnected diagrams, making the latter 1/N suppressed as compared to the fully connected diagrams. The channel where the degeneracy expected at large N is most badly broken is the pseudoscalar channel, since m η m π 0 . Therefore, we expect the Π 2+2 class of diagrams to be dominated by the π 0 and η contributions, since their contributions cancel each other to a far lesser extent than for other meson pairs such as a 2 and f 2 . For instance, the empirical ratio of the two-photon widths of a 2 and f 2 are roughly as expected if one neglects disconnected diagrams (see Table II for the source of the phenomenological values), On the other hand, using the phenomenological values [37,38] Γ η γγ = 4.35 keV, Γ π 0 γγ = 7.82 eV, (41) and the fact that Γ ∝ m 3 which does not agree well with the large-N expectation and the approximation of m strange = ∞ implicitly made here by using relations derived in N f = 2 QCD. A further channel in which the large-N expectations are not well fullfilled is the scalar sector; in particular, there does not seem to be an isovector analogue of the f 0 (600) meson. However, it turns out that the scalars make an overall small contribution to the light-by-light amplitudes.
3. The case of N f = 3 QCD Here we assume exact SU(3) flavor symmetry, m u = m d = m s . Since in the previous paragraphs we have treated the case of N f = 2 QCD, corresponding to m s = ∞, one may hope that the real world lies somewhere in between these two idealized cases.
In nature, (43) It follows directly from Q u + Q d + Q s = 0 that no diagrams with a single vector current inside a quark loop occurs in the HLbL amplitude: Π 3,1 , Π 2+1+1 and Π 1+1+1+1 do not contribute. However it is useful to keep the charges generic to derive relations from flavor symmetry and large-N arguments.
In the SU(3) symmetric theory, mesons come in octets and singlets 5 . Taking the pseudoscalar sector as an example, the transition form factor of the neutral pion is given by (Q 2 u − Q 2 d )F 8 ; neglecting again the single vector current inside a quark loop, the η meson has the transition form factor 1/3( Under the same assumption, the form factor of the η has a charge dependence given by 2/3(Q 2 Thus the pole contribution of a meson octet to the four-point function of the electromagnetic current has the following dependence on the quark charges, where Π oct is the octet contribution to Π. The corresponding expression for the singlet meson reads Matching the expressions (44) and (20), we obtain the relations where we have consistently neglected the octet contribution to Π 3+1 , Π 2+1+1 and Π 1+1+1+1 .
Under the corresponding assumption for the singlet contribution, we have Combining these equations, we finally have the following expressions for the fully connected and 5 Higher representations are allowed by symmetry, but do not seem to occur in QCD at low energies. the (2+2) disconnected class of diagrams, In the last of these equations we have set Q u = 2/3, Q d = Q s = −1/3.

C. Lattice calculation of the fully-connected vector four-point function
We now describe our method to calculate the six contractions that have fully connected quark lines (leftmost topology in Fig. 1), whereas the dominant class of disconnected diagrams (second diagram topology from the left in Fig. 1) is discussed in the next section.
We discretize the Euclidean four-point function using the local and conserved currents, as well as a contact operator, is the quark charge matrix, U µ (X) are the gauge links, and Z V is the renormalization factor for the local current. We use one local and three conserved currents. In position space, the lattice four-point function is given by The contact terms are present when two or three conserved currents coincide and serve to ensure that the conserved-current Ward identities hold. Using the backward lattice derivative ∆, these take the form In momentum space, we evaluate the Euclidean four-point function as The fully-connected contribution to Eq. (55), which is the part proportional to Tr (Q 4 ) = 17 81 , is evaluated using the method of sequential propagators. First, a point-source propagator where S is the single-flavor all-to-all quark propagator (degenerate for u and d), is computed from the origin. To concisely describe the sequential propagators, we introduce the "insertions" J µ (X) and T µ (X) for the conserved vector current and contact operator: Formally, these objects have the same size as an all-to-all quark propagator, but they are exactly zero for all sites except those that are removed from X by one lattice spacing. The point-source propagator is then combined with a plane wave and the conserved vector current insertion to form the source for new (sequential) propagators, These, in turn, are used to form sources for double-sequential propagators Finally, noting that γ 5 J µ (X) is anti-Hermitian and γ 5 T µ (X) is Hermitian, the fully-connected four-point function is obtained 6 as where . . . U denotes the expectation value over gauge fields. The sequential propagators depend on Q 1 , so a separate calculation must be done for each Q 1 . However, none of the sources for the propagators depend on Q 2 ; therefore, we are able to efficiently evaluate Π E,conn µνρσ (Q 1 , Q 2 ) for all Q 2 available on the lattice. In momentum space, the conserved-current Ward identities take the formQ whereQ µ ≡ 2 a sin aQµ 2 . We have verified that in our implementation these hold on each gauge configuration.

D. Lattice calculation of the (2+2) disconnected four-point function
We also calculate one class of disconnected diagrams to obtain an indication of their relevance. Based on the charge factor and the arguments given in section III B, the second class in Fig. 1, which we call (2 + 2) and is proportional to Tr (Q 2 ) 2 = 25 81 , is the most important. We evaluate this class of diagrams using a different lattice expression that has two local and two conserved currents: where (i, j) = (l, c) or (c, l) depending on the contraction, chosen such that each quark loop contains one local and one conserved current. Specifically, denoting the Y -to-all propagator as where is the QCD-connected expectation value over gauge fields. Each trace corresponds to a quark loop in Fig. 1.
Since X 1 and X 3 are summed over, we evaluate the traces involving S X 1 and S X 3 stochastically. To do this, we introduce a color triplet, scalar noise field φ a (X) with randomly chosen We use this as the source for two quark propagators, where each spin component is solved independently using the same noise source φ [39]. With these, we use the one-end trick [40] and obtain We reduce this stochastic noise by averaging over four noise fields per configuration, as well as using color dilution [41,42] and hierarchical probing [43] with 32 Hadamard vectors. We find that for these two-point loops hierarchical probing has no benefit over using additional noise fields; however, the noise-source propagators can be reused for one-point loops relevant for the hadronic vacuum polarization and for the other disconnected four-point diagrams, and for those loops it is beneficial. We also find that it can (if possible) in some cases be beneficial to average over the exchange of the local and conserved currents in a quark loop. We further reduce gauge noise by translating the origin of the point-source propagator S 0 and averaging over 128 point sources per gauge configuration.

A. Lattice setup
The four-point correlation functions are computed on a subset of the N f = 2 CLS (Coordinated Lattice Simulations) ensembles generated using the plaquette gauge action for gluons [45] and the O(a)-improved Wilson-Clover action for fermions [46] with the non-perturbative parameter c SW [47]. The fermionic boundary conditions are periodic in space and antiperiodic in time. We consider two different values of the lattice spacing and different pion masses in the range from 190 to 440 MeV. The parameters of the ensembles used in this work are summarized in Table I. For each ensemble, the connected four-point correlation function is computed at a few values of Q 1 = (n · 2π/T, 0, 0, 0), the first-listed component corresponding to the time direction, with Table I: Parameters of the simulations: the bare coupling β = 6/g 2 0 , the lattice resolution, the hopping parameter κ, the lattice spacing a in physical units extracted from [44], the pion mass m π , the rho mass m ρ and the number of gauge configurations. n = 1, 2, 3 on ensemble E5; n = 1, 3 on F6 and F7; n = 1, 4 on G8 and n = 2 on N6. For the (2+2) diagrams, we use n = 2 (E5) and n = 3 (F6). Then, for each value of Q 1 , the four-point correlation function is evaluated for many different values of Q 2 , corresponding to different values of Q 2 2 and ν. For the fully-connected diagrams, we used two source positions on ensembles E5, F6, F7; one source position on N6; and eight source positions on G8. In the latter case, we used the truncated-solver method [48] for the eight sources and a computation with exact inversions of the Dirac operator for bias correction on one source.
To estimate the even subtracted amplitudes, we compute the subtraction term directly at ν = 0, i.e., with Q 2 orthogonal to Q 1 . For the odd subtracted amplitudes, we use the approximation where ν 1 is the smallest available nonzero value of ν. In both cases we linearly interpolate the subtraction term in Q 2 2 to match the value in the unsubtracted term. In all tables and figures, our results for the HLbL amplitudes are multiplied by a factor of 10 6 for better readability. Since all amplitudes vanish in the limit of either Q 1 or Q 2 → 0, the signal deteriorates at small Q 2 1 (for fixed Q 2 2 ) as can be seen by comparing the left and right panels of Fig. 3.

C. Disconnected contribution to the forward light-by-light amplitudes
We now come to our results for the (2+2) disconnected diagram contribution to the eight subtracted amplitudes. We obtain this contribution with a reasonable statistical precision; however, some of the amplitudes are significantly different from zero when Q 2 2 = 0, as shown in the left panel of Fig. 2. In infinite volume, the Euclidean four-point function should vanish at this kinematic point, since a conserved current can be written as the divergence of a tensor field, J µ (x) = ∂ ν (x µ J ν (x)), so that d 4 x J µ (x) is a pure boundary term, which vanishes in the presence of a mass gap. Therefore this is a sign of significant finite-volume effects. The bulk of the effect may be removed when subtracting the amplitude at ν = 0, but some of it may remain. Figure 2 also shows that due to correlations, the subtraction significantly reduces the statistical uncertainty. The full set of subtracted amplitudes on ensemble F6 is shown in Fig. 5.   [37], as well as by [49] for the two-photon width of the f 2 (1270) meson and [38] for the π 0 width. In the case of the axial-vector mesons, the indicated width is the effective width defined in Eq. (77) and obtained phenomenologically in [6]. A cross indicates an absent or imprecise value in the PDG. An asterisk means that we use the isoscalar result divided by a factor 25/9 as explained in section III B.
A. Model description and particle content In this section, we describe how we model the hadronic γ * γ * -fusion cross section. We represent it as a sum of contributions from charge-conjugation even mesonic resonances produced in the s-channel. Specifically, we include the pseudoscalar (J PC = 0 −+ ), scalar (J PC = 0 ++ ), axial-vector (J PC = 1 ++ ) and tensor (J PC = 2 ++ ) mesons. Table II lists the most relevant light mesons with these quantum numbers. In our implementation, we limit ourselves to the lightest state in each symmetry channel. The assumption that those states are sufficient to saturate the sum rules is motivated by the fact that, at small energies, higher mass singularities are suppressed in Eq. (9). Moreover, we have revised the model used in [8] to better account for the fact that we perform fits to the fully-connected diagrams. Rather than including isovector and isoscalar mesons, we consider only isovector mesons, enhanced by a factor 34/9: we refer the reader to section III B for a justification of this approximation, which we expect to be superior. The procedure mostly modifies the contribution of the pseudoscalar sector, due to the large mass difference between the pion and the η meson. Also, since lattice simulations are performed using N f = 2 dynamical quarks, we do not include the η meson. Finally, we include the Born approximation to the γ * γ * → π + π − cross section using scalar QED, as described in Ref. [5], using a monopole vector form factor, the monopole mass being set to the ρ meson mass. Explicit formulae for cross sections used in our model are given in Appendix A. The individual contributions to the eight amplitudes from each channel are summarized in Table III.

B. Assumptions on masses and resonances
Our lattice simulations are performed at larger-than-physical quark masses. For each ensemble, the pion and ρ meson masses are determined from the pseudoscalar and vector two-point correlation functions respectively; see Table I for the obtained values. To obtain an estimate of the lowest-lying meson mass m X in every other symmetry channel, we assume that m X admits a constant additive shift relative to its physical value m phys X . The shift δm is determined from the difference between the ρ mass computed on the lattice and its experimental value, In section VI, we will test the sensitivity of our results to variations of δm by a factor of two. As for resonances, we assume that their contributions are well approximated by Breit-Wigner distributions and use the following formal substitution in the cross sections given in Appendix A, where m X and Γ X are the mass and the total width of the particle respectively. However, the remaining part of the cross section is still evaluated at s = m 2 X . For the (very narrow) pseudoscalar mesons, one can perform the integration explicitly and obtain the following contribution to the sum rules (using δ(ν − ν P ) = 2δ(s − s P ), where ν P = 1 2 (m 2 P + Q 2 1 + Q 2 2 )) : in the even case, and in the odd case, where X P ≡ ν 2 P − Q 2 1 Q 2 2 .

C. Parametrization of the form factors
In this subsection, we briefly review the available information on the transition form factors of the exchanged mesons in the hadronic model, and present the parametrization we use in fitting the lattice HLbL amplitudes. While detailed information is available in the case of the pion from lattice QCD, no experimental data is presently available at doubly virtual kinematics in any channel. In these cases, a monopole or dipole ansatz, in which the Q 2 1 and Q 2 2 dependence factorizes, is made to describe the photon-virtuality dependence, even though such an ansatz might not have the asymptotic behavior predicted by the operator-product expansion. Our motivation is that this type of parametrization is used in model calculations of a HLbL µ . Also, given our goal of performing fits to the HLbL amplitudes computed on the lattice, the number of free parameters characterizing the transition form factors should be commensurate with the precision of the lattice data.

Pseudoscalar mesons
For pseudoscalar mesons, experimental data are available when at least one photon is on-shell, and in this case a good parametrization of the data is obtained using a monopole form factor [50][51][52][53]. However, as shown in Ref. [24], a monopole form factor failed to reproduce the lattice data in the doubly-virtual case, in contrast to the LMD+V model. Furthermore, the LMD+V model is compatible with the Brodsky-Lepage behavior [54][55][56] in the singly-virtual case and with the operator-product expansion (OPE) prediction [57,58] in the Q 2 1 = Q 2 2 doubly-virtual case. We therefore use this model for the pion transition form factor, of which the parameters were determined in Ref. [24] for each ensemble listed in Table I.

Scalar mesons
Scalar mesons can be produced by two transverse (T) or two longitudinal (L) photons. Correspondingly, the amplitude is parametrized by two form factors, F T Sγ * γ * and F L Sγ * γ * . Only the first one has been measured experimentally: this was done for the f 0 (980) meson in the region Q 2 < 30 GeV 2 by the Belle Collaboration [59], and the results are compatible with a monopole form factor with a monopole mass M S = 0.800(50) GeV. Therefore, we assume the form The normalization is obtained from the experimentally measured two-photon decay width Γ γγ given by (see Table II) while the monopole mass M S will be treated as a free parameter.

Axial mesons
For axial mesons, we have two form factors, F Aγ * γ * and F Aγ * γ * , corresponding to the two helicity states of the meson. We use the same parametrization as in Ref. [5], inspired by quark models, in which 2ν = m 2 A + Q 2 1 + Q 2 2 with m A the meson mass, and assuming factorization such that A(Q 2 1 , Q 2 2 ) = A(Q 2 1 , 0)A(0, Q 2 2 )/A(0, 0) = A(Q 2 2 , Q 2 1 ). In particular, the form factor F (1) Aγ * γ * is not symmetric in the photon virtualities Q 2 1 , Q 2 2 . These form factors have been measured by the L3 Collaboration for one real and one virtual photon in the region Q 2 < 5 GeV 2 [60,61] for the isoscalar resonance. Using the previous parametrization, the authors obtain the dipole mass M A = 1040(78) MeV for the f 1 (1285) meson. We obtain the normalization of the form factors from the values given in [6] for the effective two-photon width, defined asΓ and we will consider M A as a free parameter in our fits.

Tensor mesons
We now turn our attention to the tensor mesons. The singly-virtual form factors of the isoscalar resonance f 2 for helicities Λ = 2, 1, (0, T ) have also been measured experimentally in the region Q 2 < 30 GeV 2 by the Belle Collaboration [59], where the data are compatible with a dipole form factor [6]. Therefore, we use the following parametrization for all helicities Λ = (0, T ), (0, L), 1, 2, where we allow for a different dipole mass for each helicity. The normalization of the transverse form factors is computed from the experimentally measured two-photons widths [37], Γ γγ = Γ (0) γγ + Γ (2) γγ , assuming that the ratio of helicity 2 to helicity 0 decays is r = 91.3 % (see Ref. [62]): In Ref. [6], the authors obtain the normalization of the two other form factors by saturating two different sum rules involving one real and one virtual photon; their results are summarized in Table IV.
Finally, based on large-N arguments reviewed in section III B, we assume the following relationship between the two-photon decay widths of the isoscalar and isovector mesons, In particular, we observe that this approximation works well for the tensor meson, where the two-photon decay widths have been measured both for the isovector and isoscalar mesons (see Table II).

A. Preliminary checks
In this section, we fit simultaneously the eight forward light-by-light amplitudes using the phenomenological model described in Sec. V. We have checked that we can reproduce the results given in Refs. [5,6] in the limit where only one photon is virtual to the quoted accuracy 7 (Tables I  and II of [5] and Table III  Moreover, these last three amplitudes also depend on the longitudinal scalar form factor and on the tensor form factor with helicity Λ = (0, L) which are unknown from experiment and for which we use values from phenomenology (see Table IV). As shown in the last row of Table VI, the contribution from scalar QED is always small and therefore we do not try to fit the associated monopole mass which is explicitly set to the rho mass computed on the lattice. We therefore have six fit parameters, which correspond to the monopole and dipole masses of the scalar (M S ), axial (M A ) and tensor (M for G8 are shown in appendix B, 10). The quoted error on the fit parameters is only statistical and estimated using the jackknife method. The quoted χ 2 correspond to uncorrelated fits. The χ 2 per degree of freedom are slightly above unity, with the exception of the value for ensemble E5. Here we attribute its large value to the fact that the statistical errors are smallest on E5 and that finite-volume effects could be significant for this ensemble. Given that lattice artifacts and finite-size effects are not taken into account by the χ 2 , we consider the obtained description of the data on the other ensembles to be satisfactory.
In Table VI,     Then, each row corresponds to a new fit using p ± δp and varying only one parameter at a time: the quoted number is the shift observed for the monopole/dipole mass, in units of GeV. A cross indicates that the parameter remains unchanged to all indicated digits of the central values in the first row. For instance, using Γ γγ (a 0 ) + δΓ γγ (a 0 ) instead of Γ γγ (a 0 ), the scalar monopole mass is shifted by −0.09 GeV, the other monopole/dipole masses being unaffected. In the last row, the mass shift δm applied to the spectrum (see Eq. (67)) is varied by a factor of two.  (7) 1.15  In the previous fit, only the monopole and dipole masses entering the form factors were considered as fit parameters. The other parameters (p = Γ, Γ γγ , δm, . . . ) were fixed using phenomenology as described in Sec. V. However, these parameters are sometimes associated with relatively large experimental errors (δp) or modelled (like the global mass shift in the spectrum where we assume m X = m exp X + δm with δm = m lat ρ − m exp ρ ). Therefore, we perform exactly the same fit as in the previous section but using p ± δp instead of p (and varying only one parameter at a time). In this way, we can see the influence of these parameters on the monopole and dipole masses obtained in the previous section. The results are summarized in Table VII for the ensemble F6. In this table, δm corresponds to the global mass shift applied to the spectrum (see Eq. (67)), and is multiplied or divided by a factor of two.
We observe that the experimental error on the total decay widths of the particles have a negligible effect. Increasing the two-photon width (or equivalently, the normalization of the form factor) tends to reduce the associated monopole or dipole mass. Finally, increasing the global mass shift by a factor two leads to a noticeable change in the monopole and dipole masses with little change in the χ 2 .
Varying the normalization of the form factor F (0,L) T γ * γ * leads to negligible changes in all parameters but M and M LL , which are less precisely determined on the lattice. In particular the fit is not able to determine both the dipole mass and the normalization independently, and they are highly correlated. To illustrate this point, we use the previously obtained best fit parameters and

E. Chiral extrapolations
Finally, we perform a chiral extrapolation for the six monopole and dipole masses that we have fitted to the lattice four-point function. For each of these parameters, we assume a linear dependence on m 2 π . We perform two sets of fits, either including or excluding the ensemble E5, which has the largest pion mass. The lattice results are given in Table VIII and depicted in Fig. 9 together with the fits excluding E5. The displayed errors are purely statistical. The blue points in Fig. 9 represent the ensemble N6 and therefore correspond to a finer lattice spacing than the other data points. We remind the reader that we have included only isovector mesons in the description of the fully connected amplitude, therefore all our fitted form factor parameters correspond to isovector mesons. While the results are quite stable under including or excluding ensemble E5, we consider the latter to be our final results, mainly because the χ 2 per degree of freedom of the global fit was unacceptably large for E5. We make the following observations: • The monopole mass of the scalar meson transition form factor does not depend strongly on the pion mass. After a mild extrapolation, we obtain M S = 1.04 (14) GeV at the physical pion mass. The result lies above the experimental result M S = 0.796(54) GeV from the Belle Collaboration for the isoscalar scalar meson [59].
• The axial dipole mass is also very weakly dependent on the pion mass. We obtain M A = 1.32 (7) GeV at the physical pion mass. The finer ensemble N6 suggests a value 10% larger. More ensembles would be needed to confirm whether M A is afflicted by large discretisation effects. For comparison, the L3 Collaboration obtained a dipole mass M A = 1.040(80) GeV for the isoscalar partner f 1 (1285); the measurement relied on single-virtual measurements only [60,61]. The difference in the kinematics at which the form factor was probed could be part of the reason we found a larger dipole mass, in addition to a potential genuine difference between the isospin partners. We also recall that the transition form factors have been parametrized in a fairly simplistic way (see Eq. (76) and above).
• Finally, for the tensor meson a 2 , linear extrapolations in m 2 π yield the results given in Table VIII. Fits to experimental data on the single-virtual form factor [6,59] yielded smaller values for the f 2 (1270) meson. For instance, our result for the helicity-2 transition form factor, M   , where we find that the form factors fall off far more slowly; this discrepancy could be due to the use of the factorization assumption for the dependence on the photon virtualities. On the other hand, we find agreement within the uncertainties for the scalar monopole mass and the helicity-two form factor of the tensor meson.

VII. STUDY OF DISCONNECTED DIAGRAMS
In this section, we test whether the hadronic model of section V together with the arguments summarized in section III B is consistent with the (2+2) disconnected diagrams, which we have computed on two lattice ensembles. The arguments, based on the large-N motivated idea that an isolated vector current insertion in a fermion loop gives a suppressed contribution, lead to the conclusion that the (2+2) disconnected class of diagrams contains all of the contributions from flavor-singlet meson poles, while the mesons in the adjoint representation of the flavor symmetry group contribute with a negative weight factor; the latter is (−25/9) in the SU(2) flavor case and (−2) in the SU(3) flavor case. The generic large-N expectations would further lead to the stronger conclusion that, in each J P C sector, the non-singlet resonances cancel the contribution of the flavor-singlet resonances. One channel, however, where the degeneracy is badly broken is the pseudoscalar sector, since the pion is much lighter than the η meson. Therefore, in the two-flavor theory we expect the (2+2) disconnected class of diagrams to be given to a good approximation by .. .. .
We have calculated the π 0 → γ * γ * transition form factor on the same lattice ensembles as used here in a previous publication [24]. For the η , the two-photon decay width is fairly well known experimentally; thus, assuming a vector-meson-dominance model for the virtuality dependence of the η transition form factor and using the known ρ mass on each lattice ensemble, Eq. (81) provides a prediction for the M (2+2) ..
amplitudes. We use the vector mass given in Table I. In Fig. 5 we display the prediction for the three largest subtracted amplitudes M We have obtained some evidence from our N f = 2 lattice data that the (2+2) class of diagrams is dominated by the pseudoscalar exchanges with the weight factors derived in section III B. Using the results from [67], we can now estimate the importance of the (2+2) disconnected class of diagrams in a HLbL µ in two limits: • m s = ∞, which corresponds to the two-flavor theory; • m s = m ud , which corresponds to the SU(3)-flavor symmetric theory.
We expect the real world to lie between these two predictions. In these two limits, we obtain We have used the LMD+V result for the pion (62.9 · 10 −11 ) and the VMD results for the η and η (respectively 14.5 · 10 −11 and 12.5 · 10 −11 ) quoted in [67] (see also References therein) and assigned to each contribution an uncertainty of 15%. For comparison, the N f = 2 lattice calculation [24] of the pion transition form factor and its parametrization by the LMD+V model led to the value a HLbL,π 0 µ = (65.0 ± 8.3) · 10 −11 . Taking in addition the result a HLbL µ ≈ (102 ± 39) · 10 −11 from a model calculation [68], the generic large-N based expectations imply the following estimate for the fully connected class of diagrams, These estimates give an idea of what to expect in forthcoming lattice calculations. We remark, as also pointed out in [67], that the VMD model for the η and η transition form factors is not tested in the doubly virtual case; and that the VMD form factor falls off as (Q 2 ) −2 in the limit of two large spacelike virtualities Q 2 1 = Q 2 2 = Q 2 , whereas the operator-product expansion predicts a 1/Q 2 fall-off. Thus the η and η contributions above could be somewhat underestimated due to the use of the VMD model. In the case of the pion, the 'bias' from using the VMD is −10%, relative to using the more sophisticated LMD+V model.
The only lattice calculation [19] to have presented results for a HLbL,(4) µ and a HLbL,(2+2) µ found respectively 116.0(9.6) and −62.5(8.0) in units of 10 −11 . We conclude that either these lattice results are severely underestimated, which could be due to discretization and finite-volume effects; or the hadronic model based on resonance exchanges is not viable; or the large-N inspired approximations made to estimate (82) and (83) are inadequate; or a combination of the above. A new high-statistics lattice calculation of a HLbL,(2+2) µ in a large volume would be particularly illuminating to resolve the issue, since the prediction (82) is relatively clear-cut.

VIII. CONCLUSION
With the hadronic light-by-light contribution to the muon anomalous magnetic moment a HLbL µ in mind, we have studied the eight forward light-by-light amplitudes for spacelike photons in N f = 2 lattice QCD. Via dispersive sum rules, we have tested whether the type of hadronic models used to estimate a HLbL µ provides a good description of lattice results. All in all, we found that by fitting the virtuality dependence of six meson transition form factors, we were able to describe the lattice data within statistical uncertainties. The monopole and dipole masses parametrizing the transition form factors compare reasonably well in magnitude with phenomenological determinations for the I = 0 isospin partner, with the notable exception of the dipole masses of the tensor meson for helicities Λ = 1 and Λ = (0, T ), where we find that the form factors fall off far more slowly. The simultaneous fit to all eight amplitudes allowed us to test the individual relevance of the various resonance contributions, given that they appear with different weights and signs in different amplitudes. Thus our study provides evidence, by a completely independent method, that the resonance-exchange model widely used in calculating a HLbL µ is not missing a large contribution. The (2+2) disconnected class of diagrams was computed on two lattice ensembles. We found that a parameter-free prediction based on a specific large-N argument presented in detail in section III B (see also the earlier [26]), which expresses this set of diagrams in terms of the pseudoscalar mesons alone, was compatible with the lattice data, albeit within large relative errors. Motivated by this observation, we estimated what values a lattice calculation would have to obtain for the fully connected and (2+2) set of disconnected diagrams if it is to reproduce the current model estimates of a HLbL µ . While we laid out many technical details of the method, we regard the present calculation as exploratory, and leave a more quantitative comparison of monopole and dipole masses, including an estimate of systematic errors, for the future. Indeed we were only able to perform stable fits by making model assumptions, for instance about the masses of the lightest resonances in the scalar, axial-vector and tensor sectors in N f = 2 QCD at non-physical quark masses. In addition to neglecting the three classes of diagrams containing at least one isolated vector current insertion in a quark loop, we had to assume various relations between the two-photon decay widths of isospin-partner resonances that are justified only for a large number of colors N . Also, the employed parametrization of the axial-vector resonance form factors is a further vulnerable assumption.
In the future, it would be useful to repeat the calculation of the forward light-by-light amplitudes with higher statistics, on ensembles including also the dynamical strange quark effects, and with a lighter pion mass. Especially at virtualities 0.1 GeV 2 , which can contribute significantly to a HLbL µ [67], smaller statistical errors would be beneficial to test the hadronic model more stringently. Finite-volume effects could not be addressed in any detail here, and a dedicated study would be important to carry out, given the long-range nature of the neutral pion contribution [21,22].

Notations
The metric tensor of the subspace orthogonal to q 1 and q 2 is given by such that R µν q ν i = 0 for i = 1, 2. It satisfies R µν = R νµ , R µ µ = 2 and R µ α R αν = −R µν . We use the 'mostly minus' metric convention. The virtual photon flux factor is defined through X = (q 1 ·q 2 ) 2 −q 2 1 q 2 2 = ν 2 − Q 2 1 Q 2 2 with the crossing-symmetric variable ν given by ν = q 1 · q 2 The vectors k i are defined by and satisfy k 2 i = 1, k i · q i = 0. Finally, the helicity amplitudes for the γ * (λ 1 , q 1 )γ * (λ 2 , q 2 ) → X(p X ) fusion process are related to the Feynman amplitudes by

Scalar mesons
The transition γ * (q 1 , λ 1 ) + γ * (q 2 , λ 2 ) → S where S is a scalar state can be parameterized by one transverse (F T Aγ * γ * ) and one longitudinal (F L Aγ * γ * ) form factor and is described by the following matrix element The only non-zero helicity amplitudes are given by The two-photon decay width is given by and from Eqs. (2) and (8) (A10)