Hadronic light-by-light contribution to (g−2 ) μ from lattice QCD with SU(3) flavor symmetry
En-Hung Chao, Antoine Gérardin, Jeremy R. Green, Renwick J. Hudspith, Harvey B. Meyer
MMITP/20-036CERN-TH-2020-109
Hadronic light-by-light contribution to ( g − µ from lattice QCD with SU(3) flavor symmetry En-Hung Chao, Antoine G´erardin, Jeremy R. Green, Renwick J. Hudspith, and Harvey B. Meyer
1, 41
PRISMA + Cluster of Excellence & Institut f¨ur Kernphysik,Johannes Gutenberg-Universit¨at Mainz, D-55099 Mainz, Germany Aix Marseille Univ, Universit´e de Toulon, CNRS, CPT, Marseille, France. Theoretical Physics Department, CERN, 1211 Geneva 23, Switzerland Helmholtz Institut Mainz, Johannes Gutenberg-Universit¨at Mainz, D-55099 Mainz, Germany (Dated: July 15, 2020)We perform a lattice QCD calculation of the hadronic light-by-light contributionto ( g − µ at the SU(3) flavor-symmetric point m π = m K (cid:39)
420 MeV. The rep-resentation used is based on coordinate-space perturbation theory, with all QEDelements of the relevant Feynman diagrams implemented in continuum, infinite Eu-clidean space. As a consequence, the effect of using finite lattices to evaluate theQCD four-point function of the electromagnetic current is exponentially suppressed.Thanks to the SU(3)-flavor symmetry, only two topologies of diagrams contribute,the fully connected and the leading disconnected. We show the equivalence in thecontinuum limit of two methods of computing the connected contribution, and intro-duce a sparse-grid technique for computing the disconnected contribution. Thanksto our previous calculation of the pion transition form factor, we are able to correctfor the residual finite-size effects and extend the tail of the integrand. We test ourunderstanding of finite-size effects by using gauge ensembles differing only by theirvolume. After a continuum extrapolation based on four lattice spacings, we obtain a hlbl µ = (65 . ± . ± . × − , where the first error results from the uncertain-ties on the individual gauge ensembles and the second is the systematic error of thecontinuum extrapolation. Finally, we estimate how this value will change as thelight-quark masses are lowered to their physical values. a r X i v : . [ h e p - l a t ] J u l I. INTRODUCTION
Electrons and muons carry a magnetic moment aligned with their spin. The proportion-ality factor between the two axial vectors is parameterized by the gyromagnetic ratio g . InDirac’s theory, g = 2, and for a lepton family (cid:96) one characterizes the deviation of g from thisreference value by a (cid:96) = ( g − (cid:96) /
2. Historically, the ability of Quantum Electrodynamics(QED) to quantitatively predict this observable played a crucial role in establishing quantumfield theory as the framework in which particle physics theories are formulated.Presently, the achieved experimental precision on the measurement of the anomalousmagnetic moment of the muon [1], a µ , is 540 ppb. At this level of precision, such a measure-ment tests not only QED, but also the effects of the weak and the strong interaction of theStandard Model (SM) of particle physics. Currently there exists a tension of about 3.7 stan-dard deviations between the SM prediction and the experimental measurement. The statusof this test of the SM is reviewed in [2–5]. At the time of writing, the E989 experiment atFermilab is performing a new direct measurement of a µ [6], and a further experiment usinga different experimental technique is planned at J-PARC [7]. The final goal of these experi-ments is to reduce the uncertainty on a µ by a factor of four. A commensurate reduction ofthe theory error is thus of paramount importance.The precision of the SM prediction for a µ is completely dominated by hadronic uncertain-ties. The leading hadronic contribution enters at second order in the fine-structure constant α via the vacuum polarization and must be determined at the few-permille level in orderto match the upcoming precision of the direct measurements of a µ . The most accuratedetermination comes from the use of e + e − → hadrons data via a dispersion relation, al-though lattice QCD calculations have made significant progress in computing this quantityfrom first principles [5, 8]. The hadronic light-by-light (HLbL) scattering contribution a hlbl µ ,which is of third order in α , currently contributes at a comparable level to the theory un-certainty budget and is being addressed both by dispersive and lattice methods; see [9–11]and references therein.Our approach for determining a hlbl µ is based on coordinate-space perturbation theorywhere the QED elements of the Feynman diagrams yielding a hlbl µ are precomputed in infi-nite volume, and only the four-point amplitude of the electromagnetic current is actuallycomputed on the lattice. Here we compute the lattice contribution at a point in the spaceof light quark masses corresponding to QCD with exact SU(3)-flavor (denoted SU(3) f inthe rest of this work) symmetry. Furthermore, the sum of the three light quark masses isapproximately the same as in nature. These two conditions leads to a degenerate mass ofpions, kaons and the eta meson of about 420 MeV.Our motivation for calculating a hlbl µ at the SU(3) f -symmetric point is twofold. First, thelattice calculation itself is simplified in that only two out of five classes of Wick contractionscontribute, due to the vanishing trace of the quark electric charge matrix. In addition, theoverall lattice calculation is computationally far cheaper than for physical quark masses, sothat more tests of systematic errors can be performed. Second, the interpretation of theresults is simplified: the SU(3) f -symmetry reduces the number of unknown parameters inmodel estimates based on the exchange of the lightest mesons. In particular, the transitionform factor (TFF) of the pion, which describes the coupling of the neutral pion to twovirtual photons, has been calculated [12] on the lattice ensembles that we use. The TFFof the eta meson coincides with the TFF of the pion up to a simple overall charge factor.Of the pseudoscalar mesons, only the TFF of the η (cid:48) remains independent and is largely x y z yz x y zx y zx yzx FIG. 1. Five Wick-contraction topologies that are necessary for the calculation of a hlbl µ as listedin Table I. For each diagram, nontrivial permutations of the four quark-photon vertices yield thenumber of contractions listed in that table. The two diagrams in the first line are the dominantones and the other three vanish in the SU(3) f limit. unknown at the SU(3) f -symmetric point, however experimental information is available forthe two-photon decay width (which provides the coupling strength to two real photons) andsome experimental results are available for the singly as well as the doubly-virtual formfactor [13–15], although only for relatively large virtualities above 1.5 GeV .The simplified connection to model estimates enables our work to provide a valuablecross-check for the predictions of hadronic models and dispersive methods; this work isthereby complementary to lattice calculations directly aiming at a hlbl µ for physical quarkmasses [16]. At the same time, this study allows us to learn about the size of various sourcesof systematic error, particularly the finite-size effects, and how well we are able to correctfor them semi-analytically.The rest of this paper is organized as follows: We begin by presenting our methodology inSection II, including the two methods we will investigate for computing the quark-connectedcontribution to a hlbl µ . Section III begins with a description of the lattice ensembles used inthis work, as well as an example of the lowest-lying relevant meson spectrum for one of ourensembles. We then discuss the lattice determination of the π and η transition form factorsused for the modelling and finite-size correction of our data. Section IV discusses the variousmodel predictions for the integrand at the SU(3) f -symmetric point and confronts these witha selection of our lattice data. Results for the fully-connected class of Wick contractions arepresented in Section V, and those for the non-vanishing disconnected class in Section VI. Inboth cases, the lattice results are compared to the prediction for the exchange of pseudoscalarmesons at the integrand level. The main result of this paper — a hlbl µ at the SU(3) f -symmetricpoint, Eq. (31) — is obtained in Section VII, which also contains a discussion of how thisresult will change for physical values of the quark masses. We summarize our findings andconclude in Section VIII, and various technical aspects of the calculation are described inmore detail in the appendices (A, B, C). II. METHODOLOGY
One can compute the light-by-light scattering contribution to the g − a hlbl µ = m µ e (cid:90) d y (cid:90) d x ¯ L [ ρ,σ ]; µνλ ( x, y ) i (cid:98) Π ρ ; µνλσ ( x, y ) , (1)where ¯ L is a QED kernel and i (cid:98) Π is a spatial moment of the connected Euclidean four-pointfunction in QCD, i (cid:98) Π ρ ; µνλσ ( x, y ) = − (cid:90) d z z ρ (cid:101) Π µνσλ ( x, y, z ) , (2) (cid:101) Π µνσλ ( x, y, z ) ≡ (cid:68) j µ ( x ) j ν ( y ) j σ ( z ) j λ (0) (cid:69) QCD , (3)with j µ ( x ) the hadronic component of the electromagnetic current j µ ( x ) = 23 ( uγ µ u )( x ) −
13 ( dγ µ d )( x ) −
13 ( sγ µ s )( x ) . (4)The QCD correlation function consists of all the various ways one can contract four vectorcurrents, as shown in Fig. 1, all of which are “disconnected” except for the fully-connectedcontribution. At the flavor-symmetric point, only the upper two topologies contribute. Awayfrom this point it is expected, from large- N c arguments as well as from numerical evidenceby the RBC/UKQCD collaboration [9], that the remaining topologies are suppressed. Thenumber of contractions for each topology is given in Table I. conn 2 + 2 3 + 1 2 + 1 + 1 1 + 1 + 1 + 1 m l (cid:54) = m s m l = m s More information on the infinite-volume QED kernel ¯ L ( x, y ) can be found in [17]. Wemake use of O (4) symmetry to simplify the integral further, a hlbl µ = (cid:90) ∞ d | y | f ( | y | ) , (5)where in the starting representation (1), f ( | y | ) = m µ e π | y | (cid:90) d x ¯ L [ ρ,σ ]; µνλ ( x, y ) i (cid:98) Π ρ ; µνλσ ( x, y ) . (6)We will often display this integrand and always denote it by f ( | y | ), even though in practicewe employ modified representations of a hlbl µ . One type of modification concerns the kernel¯ L [ ρ,σ ]; µνλ ( x, y ). Various subtraction terms to the kernel have been proposed to beneficiallychange the shape of the integrand [9, 18] without changing the resulting integral. Theimportance of performing such subtractions cannot be understated as the unsubtractedkernel is poorly suited for practical lattice simulations due to being too peaked at shortdistances.Here we make extensive use of a new subtraction scheme for the QED kernel [18],¯ L (Λ)[ ρ,σ ]; µνλ ( x, y ) = ¯ L [ ρ,σ ]; µνλ ( x, y ) − ∂ ( x ) µ ( x α e − Λ m µ x / ) ¯ L [ ρ,σ ]; ανλ (0 , y ) − ∂ ( y ) ν ( y α e − Λ m µ y / ) ¯ L [ ρ,σ ]; µαλ ( x, , (7)where Λ is an arbitrary, tuneable, dimensionless free parameter. When Λ = 0 we havethe ¯ L (2) kernel of [17] and as Λ → ∞ we recover the unsubtracted ¯ L (0) . The benefit ofsuch a choice of kernel is that we are able to tune the shape of the integrand to reducethe long-distance effects while still preserving the beneficial properties of short-distancesubtractions. An investigation of this kernel with the infinite-volume lepton loop is presentedin Appendix B.In this section we will write the formulae for obtaining a hlbl µ in terms of continuous in-tegrals. The lattice, however, is discrete so we can only approximate these integrals withfinite sums , (cid:90) d x ≈ a N t / − (cid:88) t/a = − N t / N z / − (cid:88) z/a = − N z / N y / − (cid:88) y/a = − N y / N x / − (cid:88) x/a = − N x / , (8)where a is the lattice spacing and N µ = L µ /a . In Eq. (2), we set z ρ = L/ → f ( | y | ) with respect to | y | can be performed with the trapezoid rule and in practice we will average over equivalentvalues of f ( | y | ) to both increase statistical precision and reduce computational cost. Lateron in the analysis we will show results for the partially-integrated value of a µ , a µ ( | y | max ) = (cid:90) | y | max d | y | f ( | y | ) (9)with the hope that we will see a plateau at large-enough | y | max , indicating that our integralhas saturated.In the following three subsections we will discuss two methods to compute the connectedcontribution: the first, Method 1, is a direct calculation of the three pairs of connected Wickcontractions; the second, Method 2, uses rearrangements of the integrand expressed in termsof only one pair of Wick contraction to make the calculation cheaper. We will then give themethodology used for the calculation of the disconnected contribution, which also uses somerearrangements in the integrand. A. Connected contribution (Method 1)
The Wick contractions for the connected contribution are shown in Fig. 2. With fourlocal vector currents, the four-point correlation function can be written as (cid:101) Π conn µνσλ ( x, y, z ) = 1881 Z V (cid:16) (cid:101) Π (1) µνσλ ( x, y, z ) + (cid:101) Π (2) µνσλ ( x, y, z ) + (cid:101) Π (3) µνσλ ( x, y, z ) (cid:17) , (10) If using open boundary conditions (as we mostly do in this work) some care in truncating the temporalsum is required as to not incorporate boundary effects. z x y z x y x y z FIG. 2. Wick contractions for the connected contribution. Each diagram represents two contrac-tions with quark flow in opposite directions. where Z V is the renormalization factor of the vector current and each term represents onediagram in Fig. 2, i.e. a pair of Wick contractions. We will use the Z V values from [19]without the O ( a )-improvement terms proportional to quark mass combinations. The value18 /
81 is the necessary charge factor for the degenerate light and strange quarks. Writingthe (cid:101) Π (1 , , s in terms of propagators yields (cid:101) Π (1) µνσλ ( x, y, z ) = − (cid:104) Tr [ S (0 , x ) γ µ S ( x, y ) γ ν S ( y, z ) γ σ S ( z, γ λ ] (cid:105) U , (cid:101) Π (2) µνσλ ( x, y, z ) = − (cid:104) Tr [ S (0 , y ) γ ν S ( y, z ) γ σ S ( z, x ) γ µ S ( x, γ λ ] (cid:105) U , (cid:101) Π (3) µνσλ ( x, y, z ) = − (cid:104) Tr [ S (0 , y ) γ ν S ( y, x ) γ µ S ( x, z ) γ σ S ( z, γ λ ] (cid:105) U , (11)where S ( x, y ) is the quark propagator with sink at x and source at y and (cid:104)· · · (cid:105) U is theexpectation value over gauge configurations. The trace is over both Dirac and color indices.In Method 1, the six contractions are computed explicitly and we choose the direction y/a ∝ (1 , , , . In this method, we first compute the point-to-all propagator S ( · ,
0) and the sixsequential propagators using the fields z [ ρ γ σ ] S ( z,
0) as sources (the anti-symmetrization ofthe indices ρ and σ is imposed by the symmetries of the QED kernel). Then, for each valueof | y | used to sample the integrand in Eq. (5), we compute the point-to-all propagator S ( · , y )and the six sequential propagators using the fields z [ ρ γ σ ] S ( z, y ) as sources. Therefore, for N evaluations of the integrand, we need 7( N + 1) propagator inversions. In our set-up allcurrents are local except for the one at x , which will be either local or the point-split, thelatter implying a suitable modification of Eq. (11).In practice, since we are (mostly) using open-boundary conditions in the time direction,the origin is located somewhere near the middle time-slice and randomly distributed inthe spatial volume. We also average over the 16 combinations y/a = ( ± n, ± n, ± n, ± n ) toincrease statistics. B. Connected contribution (Method 2)
The idea for Method 2 is simple: we pick a reference diagram that is easiest to computeand use a change of variables in the integrals to relate the other diagrams to this reference.Here, we pick the diagram that does not have a propagator from the origin to y , or from z to x (the leftmost of Fig. 2). Such a choice avoids the extra inversions required for the Note that the fourth power in Z V in Eq. (10) comes from the use of four local currents. A conservedcurrent does not require such a factor. This should thus be adapted correspondingly in some particularcases where the effect of conserved currents is studied in this work. We use the notation ( x , y , z , t ) where time is the last component sequential sources as we discussed in Method 1. Simplistically, for two samples of a single | y | (i.e. + y and − y ) we only need two point-to-all propagators. In fact, we can do muchbetter than this if we keep propagators in memory so that they can be reused.The diagrams (cid:101) Π (2) and (cid:101) Π (3) are related to (cid:101) Π (1) by applying symmetries to Eq. (11):Euclidean-space translations and inversions, along with cyclicity of the trace [17, 18]. Thisleads to our master formula for Method 2: a conn µ = − Z m µ e π (cid:90) d | y || y | (cid:90) d x (cid:18) ( ¯ L (Λ)[ ρ,σ ]; µνλ ( x, y ) + ¯ L (Λ)[ ρ,σ ]; νµλ ( y, x ) − ¯ L (Λ)[ ρ,σ ]; λνµ ( x, x − y )) (cid:90) d z z ρ (cid:101) Π (1) µνσλ ( x, y, z )+ ¯ L (Λ)[ ρ,σ ]; λνµ ( x, x − y ) x ρ (cid:90) d z (cid:101) Π (1) µνσλ ( x, y, z ) (cid:19) . (12)In infinite volume, the integral is equal to the result of Method 1, but the integrand f ( | y | )will generally be different. As a result, the systematic effects in Method 2 will be differentfrom Method 1.We invert point-source propagators along the line ( n, n, n, n + t min /a ), and what we call | y | will be the difference between these points; here the integer n ranges from 0 to N i / t min ) is chosen to be suitably away from our (usually)open temporal boundary, ideally m π t min > m π ( L t − t max ) >
4. This line of propagatorsources was chosen as we typically have a large anisotropy in the temporal direction, allowingus to achieve sizeable values of | y | .In our implementation, we keep all N i / y and origin, so for N propagator inversions we have N ( N − | y | . To further boost statistics wewill also average other directions that give the same | y | , e.g. y/a = ( ± n, ± n, ± n, t max /a − n ). C. Disconnected contribution
The disconnected contribution can be computed from the two-point contractionΠ µν ( x, y ) = − Re (Tr[ S ( y, x ) γ µ S ( x, y ) γ ν ]) . (13)An important point to note is that Π µν has a vacuum expectation value (VeV) that must besubtracted to ensure that the two “disconnected” quark loops are still connected by gluons.To this end, we define ˆΠ µν ( x, y ) = Π µν ( x, y ) − (cid:104) Π µν ( x, y ) (cid:105) U , (14)and use this to compute the 2 + 2 quark-disconnected contribution to a hlbl µ , a disc µ = − Z m µ e π (cid:90) d | y || y | (cid:90) d x (cid:28) ( ¯ L (Λ)[ ρ,σ ]; µνλ ( x, y ) + ¯ L (Λ)[ ρ,σ ]; νµλ ( y, x )) ˆΠ µλ ( x, (cid:90) d z z ρ ˆΠ σν ( z, y )+ ¯ L (Λ)[ ρ,σ ]; µνλ ( x, y ) ˆΠ µν ( x, y ) (cid:90) d z z ρ ˆΠ σλ ( z, (cid:29) U . (15)In our implementation we compute a grid of point-source propagators from some t min to t max wrapping completely around the spatial directions alternating between (2 n, n, n, m + t min /a ) and (2 n + 1 , n + 1 , n + 1 , m + 3 + t min /a ) where n and m take values between 0and N i / −
1, and 0 and ( t max − t min ) / (6 a ) respectively. (0,0,0,0) (2,2,2,0)(1,1,1,3) (3,3,3,3)(0,0,0,6) (2,2,2,6)(1,1,1,9) (3,3,3,9) FIG. 3. Illustration of the grid of propagators used in the toy example given in the text. Filledcircles indicate sites that lie on the lattice, gray circles are ones accessible through periodicity andthe dashed lines indicate lines along our preferred directions (1 , , ,
3) and (2 , , , Giving a toy example; for a 4 ×
12 periodic lattice we would invert point-source prop-agators at (0 , , , , (2 , , , , (1 , , , , (3 , , , , (0 , , , , (2 , , , , (1 , , , , and(3 , , , y lies in a direction we like, e.g. along the directions (1 , , ,
3) and(2 , , , z and x integrations are saved to disk, the VeV is subtracted,and the final contraction of indices is performed offline in the analysis.The benefit of such a method over calculating each separate | y | individually is thequadratic growth of equivalent combinations of | y | that are available as the number ofsource fields grows. For our 4 ×
12 example we have 6 values of equivalent | y | per sourceposition. It would be impractical to perform this calculation without exploiting this fact.We find it is beneficial to truncate the x -integral in our set-up to avoid negligible, noisycontributions at large distances. As our kernels are computed on-the-fly this is also useful asa computer-time saving measure as it reduces the number of QED kernel calls. We truncatethe integral over x to points that lie within a maximum distance [( r/a ) = 81 , , , y , while we do nottruncate the z -integral. Although these truncations differ in physical volume, the ensemblesH101 and H200 were tested for smaller and larger values of ( r/a ) on a subset of our dataand even smaller values of ( r/a ) than used here, for example 49 and 64 for the ensembleH101, were found to be consistent.To reiterate, in our full calculation we do not perform the integrals for all possible y -vectors, as there are many that we expect to have bad finite-volume or discretization effects. Label β κ size bdy. cond. a (fm) m π,K,η (MeV) m π L U103 3.40 0.13675962 24 ×
128 open 0.08636(98)(40) 415(5) 4.4H101 3.40 0.13675962 32 ×
96 open 0.08636(98)(40) 418(5) 5.9B450 3.46 0.13689 32 ×
64 periodic 0.07634(92)(31) 417(5) 5.2H200 3.55 0.137 32 ×
96 open 0.06426(74)(17) 421(5) 4.4N202 3.55 0.137 48 ×
128 open 0.06426(74)(17) 412(5) 6.4N300 3.70 0.137 48 ×
128 open 0.04981(56)(10) 421(5) 5.1TABLE II. Lattice ensembles with SU(3) f -symmetry used in this work. Each ensemble isparametrized by the gauge coupling parameter β , the quark hopping parameter κ , the latticesize, and the temporal boundary condition. The lattice spacing a was determined in Ref. [21]. . . . . a [fm ]2 . . . L [f m ] U103H101B450H200N202N300
FIG. 4. Spatial extent L and lattice spacing a for the ensembles with SU(3) f -symmetry used inthis work. For instance, in the toy example we could calculate f ( | y | ) for the y -vector y/a = (0 , , , ≈
80% of thepossible y -vectors and keep only (the modulus of) those parallel to (1 , , , , , , , , , III. LATTICE PARAMETERS AND PROPERTIES OF SU(3) f -SYMMETRICQCD Calculations have been performed on lattice ensembles provided by the CoordinatedLattice Simulations (CLS) initiative [20], which have been generated using three flavorsof non-perturbatively O ( a )-improved Wilson-clover fermions and with the tree-level O ( a )-improved L¨uscher-Weisz (Symanzik) gauge action. In particular, we consider only those withSU(3) f -symmetry. On these ensembles, the mass of the octet of light pseudoscalar mesons isapproximately 420 MeV. These ensembles are summarized in Table II and in Fig. 4: thereare four lattice spacings, as well as two pairs of ensembles that differ only by their volume.0 J P C octet singlet0 − + π η (cid:48) ++ a f ++ a f −− ρ ω + − b h A. The SU(3) f meson spectrum At long distances, the leading contributions to the four-point correlator emanate fromthe lowest-lying mesons. Clearly, the degenerate π and η -pole contributions dominate atthe longest distances. However, at intermediate distances the contributions from heavierstates, including resonances in the two-pseudoscalar channel, may also play a significantrˆole. Especially since in our calculation at the SU(3) f -point the pseudoscalar mesons arenot appreciably light compared to the other mesons. Therefore, we present results from alimited spectroscopic study of the low-lying meson spectrum in all relevant channels withtotal angular momentum J ≤ C -even) mesons can couple to two electromagnetic cur-rents and thus contribute in exchange diagrams to the four-point function. However, fromthe Wick-contraction structure, one can view Method 2 as using charged vector currents(see Appendix A) and it is possible to couple a C -odd state to two of them. Therefore, theintegrand for Method 2 can receive contributions from C -odd mesons such as the rho [18].For that reason, we also inspect the spectrum of C -odd mesons, even if their contributionsto a hlbl µ must vanish in the infinite-volume limit once all integrals are performed.Our dedicated spectroscopy calculation is performed on the ensemble H101. This madeuse of the distillation framework [22] (for the connected hadron two-point functions) andits stochastic formulation [23], which was essential for efficiently computing the discon-nected hadron two-point functions needed in the singlet sector. However, we have only usedsmeared quark bilinears as interpolating operators, and the lack of non-local multimeson-like interpolators means that this calculation should not generally be considered as a robustdetermination of the spectrum. This analysis precludes the use of finite-volume quantizationconditions; therefore, only the approximate locations of resonances can be found, providedthat they are narrow. Nevertheless, since we compute diagonal correlators, the effectivemasses taken at any Euclidean time provide an upper bound on the ground-state energy ina given channel. This observation is particularly useful in the flavor-singlet 0 ++ channel,which turns out to admit a stable ground-state meson. Such a stable scalar meson hasbeen found previously [24] in a lattice calculation at a similar pion mass, though not at anSU(3) f -symmetric point.Results are shown in Fig. 5 and summarized in Table III. Note that since the signal inthe flavor-singlet sector is much worse, the plateau fits had to be done at relatively shorttime separations, so that the shown uncertainty on the mass may be an underestimate.In addition, the plateau for the flavor-octet scalar is relatively poor, which might be due1 t / a . . . . . . a m e ff octet + − ( b )1 ++ ( a )0 ++ ( a )1 −− ( ρ )0 − + ( π ) t / a . . . . . . a m e ff singlet + − ( h )1 ++ ( f )0 ++ ( f )1 −− ( ω )0 − + ( η ) − + ++ ++ −− + − . . . . . . a m octetsinglet FIG. 5. The meson spectrum of ensemble H101. Effective masses are shown for each channel inthe flavor octet (left) and singlet (center) sectors, along with bands indicating the plateau fits.Note that in many cases (in particular the octet pseudoscalar and vector channels), the error bar issmaller than the plotted symbol. Fitted masses are shown in the right panel; the outer error barsinclude an estimate of systematic uncertainty obtained by shifting the plateau fit range. Horizontaldashed lines indicate thresholds at twice and three times the octet pseudoscalar mass. to its coupling to two octet pseudoscalars in S-wave and the plateau’s proximity to thecorresponding threshold. In the case of a mis-identified plateau, the true ground-stateenergy would be lower, since the effective masses shown here are expected to approach theground-state energy from above.For the flavor-octet sector, after the π , the next-longest-distance meson-exchange con-tributions in the integrand for a hlbl µ come from the a and (for the Method 2 integrand)the ρ , both of which sit near 2 m π . For the flavor-singlet sector, the three correspondingmesons are also the lightest, although the ordering is different, with the f (or σ ) being abound state with mass near 680 MeV and the η (cid:48) mass sitting higher, somewhere near 2 m π .For the disconnected diagrams, which receive the difference between contributions from theexchanges of singlet and octet mesons (see section IV A), the pseudoscalars will providethe longest-distance contribution but we should also expect a significant contribution fromscalars. The difference between the octet and singlet vector meson masses is very small; ifthe same holds for their form factors, then their combined contribution to the integrand willbe negligible. B. Pion-pole and η -pole contributions to a hlbl µ In Ref. [12] a model independent parametrization of the pion TFF has been obtained onthe same set of lattice ensembles as used in this work. These results can be used to computethe pion-pole contribution, a hlbl; π µ , at the SU(3) f -symmetric point in the continuum limit.This result will be used in Section VII C to obtain a rough estimate of a hlbl µ at the physicalpion mass. More importantly, in Sections V A, V B, and VI, a vector meson dominance(VMD) parametrization of the TFF will be used to estimate the finite-size corrections toour lattice data at the symmetric point. For this purpose, we have performed a new fit to2 Label F π γγ (0 ,
0) [GeV − ] M V [MeV]H101 0.237(5) 921(13)B450 0.230(5) 942(25)H200 0.219(6) 979(26)N202 0.227(5) 952(15)N300 0.214(5) 1001(23)TABLE IV. The VMD fit parameters of the π transition form factor. the data from [12] assuming a VMD parametrization, where we have restricted the fit tothe singly-virtual case where the model provides a good description of our data. The fitparameters used in this work are collected in Table IV.In [12], two strategies have been used to extract the pion TFF, and the correspondingresults for the pion-pole contribution are displayed in Fig. 6 for two different discretizationsof the 3-point correlation function (see [12] for details). In the first strategy, the TFF wascomputed on each lattice ensemble separately. This allowed us to determine the pion-polecontribution for different values of the pion mass and lattice spacing. The physical value wasthen obtained by a combined chiral-continuum extrapolation of a hlbl; π µ . We have repeatedthis analysis but now restrict the fit to only the ensembles included in this study. Here weuse the pion mass of the given ensemble (instead of the physical one as done in [12]) in theweight functions that appear in Eq. (74) of [12]. This leads to the result (1) of Fig. 6.In the second strategy, the pion TFF was directly extrapolated to the physical pointusing a global fit that includes several ensembles including, and away from, the SU(3) f -point. From the resulting fit parameters, we can extract the pion TFF in the continuumlimit at a pion mass of 420 MeV. Using this result in Eq. (74) of [12], we obtain the greypoint of Fig. 6, which we find to be in very good agreement with the first estimate.The second strategy has the advantage of using an expanded set of ensembles (15 in total)to determine the TFF, a hlbl; π µ = 21 . . × − (SU(3) f point) (16)at the SU(3) f -point which is to be compared to the physical-pion value a hlbl; π µ = 59 . . × − (physical point) . (17)Unsurprisingly, we observe a strong dependence on the pion mass. The smaller pion con-tribution in Eq. (16) compared with (17) is due roughly in equal parts to the heavier pionmass and to the reduced coupling to photons, as can be seen by comparing the entries inTable IV to the physical value, F π γγ (0 , ≈ .
274 GeV − .At the SU(3) f -symmetric point, the η is mass-degenerate with the pion, m η (cid:39)
420 MeV,and contributes exactly 1 / π exchange to a hlbl µ , i.e. a hlbl; ηµ = (7 . ± . × − (SU(3) f -point) . (18)A lattice calculation at the physical point for this contribution is not yet available, but themost recent estimate , which comes from Canterbury approximants [5, 26], is a hlbl ,η, phys µ = An older estimate is a hlbl ,η, phys µ = 14 . × − using the VMD model [25]. . . . . (1) a hlbl; π µ = 19 . . × − (2) a hlbl; π µ = 21 . . × − a h l b l; π µ × a [fm ] FIG. 6. The pion-pole contribution to a hlbl µ at the SU(3) f symmetric point with a pion mass of420 MeV. Blue and red points correspond to two different discretizations of the 3-pt correlationfunction. The results (1) and (2) are explained in the main text. . × − . This comparison shows that at the physical point the η gives a much largercontribution to a hlbl µ than at the SU(3) f -symmetric point, in spite of being heavier. This canbe traced back to its much larger coupling to two photons, F ηγγ (0 , − ] (cid:39) (cid:26) .
12 SU(3) f point , .
27 physical point . (19) IV. HADRONIC MODELS VS. THE LATTICE INTEGRAND f ( | y | ) In subsection IV A, we state the theoretical predictions for the integrand f ( | y | ) corre-sponding to the quark-connected diagrams in Method 1 and Method 2, as well as for thedisconnected diagram. Then, in subsection IV B, we compare the lattice | y | -integrand ob-tained by Method 1 for the quark-connected contribution to hadronic model predictions.For this purpose, we will focus on the ensembles N202 and H200. N202 has the largestphysical volume ( L = 3 .
08 fm) of all ensembles considered in this work, and H200 (with L = 2 .
06 fm) differs only by its comparatively smaller volume. Since these only differ bytheir volume, they allow us to test our understanding of finite-volume effects. Finally, a com-parison of the theory predictions for the integrand of the quark-disconnected contributionsto the corresponding lattice data is made in subsection IV C.
A. Predictions for the integrand
In order to gain some insight into the various contributions to the quantity a hlbl µ , wewill compare predictions for the pseudoscalar exchanges as well as the contribution froma constituent-quark loop, and a charged-pseudoscalar loop, to the lattice data. Here, we4provide some details as to how these predictions are obtained. Expressions for the ampli-tudes i (cid:98) Π ρ ; µνλσ calculated with a fermion loop, a charged-pion loop, or with the pseudoscalarexchange can be found in [17, 27]. In QCD with exact SU(3) f -symmetry, the η contributionto a hlbl µ is always one third of the contribution of the π .The flavor structure of single-meson exchanges in the different Wick-contraction topolo-gies of four-point functions was discussed in detail in Ref. [28], including the case of N f = 3QCD. The quark-connected diagrams receive contributions only from flavor-octet mesons,enhanced by a factor of three. This is compensated by the quark-disconnected diagrams,which contain the differences between the flavor-singlet and twice the flavor-octet mesoncontributions. In the large- N c limit, the singlet and octet contributions to the disconnecteddiagrams will cancel as their spectra becomes the same. For QCD, this degeneracy is moststrongly broken in the pseudoscalar sector, where the octet’s mass is far below the singlet’s.To interpret the integrand in Method 2, one needs a mapping of individual quark-levelWick contractions onto meson-exchange diagrams [18] (see Appendix A for a derivationbased on partially quenched chiral perturbation theory). It turns out that (cid:101) Π (1) µνσλ ( x, y, z )does not contain the meson-exchange diagram in which the π and η propagate between thepair of vertices (0 , y ) and ( x, z ). Also, the normalization of the two other ( π , η )-exchangediagrams is such that (cid:101) Π conn µνσλ contains the same ( π , η ) contribution as (cid:101) Π µνσλ , enhanced (inthe present SU(3) f case) by the charge factor three.In addition, we need the mapping of individual quark-level disconnected diagrams ontomeson-exchange diagrams. Here it turns out that there is a one-to-one mapping between agiven quark-level diagram and a meson-exchange diagram. For instance, take a quark-leveldiagram, consisting of two quark loops each containing two vectorial vertices; each quarkloop thus defines a pair of vertices. Such a diagram is in one-to-one correspondence withthe diagram in which the meson is exchanged between the two pairs of vertices. For octetmesons, the latter diagram has a weight of − a hlbl µ is sometimes modelled by a constituent-quarkloop. This corresponds to an effective degree of freedom, and the mass assigned to the‘constituent quark’ is typically on the order of 300 MeV [31]. In the following sub-section,we will address the question “to what extent does such a contribution describe the short-distance part of the | y | integrand?”. In this case, the Wick contractions and weight factorsof the constituent quark simply correspond to those of the fundamental quark degrees offreedom.We have computed the charged pion and kaon loop contributions in the framework ofscalar QED [17, 18]. The contribution of these loops to the set of quark-connected contrac-tions and their contribution to the set of quark-disconnected contractions add up coherently,the latter being twice as large as the former. This result can again be derived in partially-quenched perturbation theory. While the ( π , η ) exchange contribution to the full a hlbl µ isthree times smaller than its contribution to a conn µ , the charged-pseudoscalar loop contributionis three times larger. The charged-pseudoscalar loop’s contribution might seem negligiblein comparison to the integrand of the quark-connected, but it need not be negligible in thefull integrand. The enhancement of the π -exchange contribution in quark-connected diagrams was first pointed outin Refs. [29, 30]. -20 0 20 40 60 80 100 120 0 0.5 1 1.5 2 f( | y | ) x [f m - ] |y| [fm] kernel L ( Λ ) , method 1, Λ =0.16 Lattice: N202 pt-splitQuark loop m q =350MeV π o + η (in finite-vol)( π ,K) loopTotal FIG. 7. Integrand for the connected contribution on ensemble N202 using Method 1 with Λ = 0 . x . The integrand is compared to the prediction forthe exchange of the π and η mesons with a VMD transition form factor, which is expected toprovide a good approximation to the tail. In addition, an attempt to describe the short-distancecontribution with a constituent-quark loop with a quark mass of 350 MeV is made. B. Lattice connected contribution
Fig. 7 displays the integrand obtained with Method 1 and Λ = 0 .
16. It is comparedto the integrand for the exchange of the ( π , η ) mesons with a VMD TFF and using theparameters of Table IV. Beyond 1.5 fm, the prediction is consistent with the lattice data,albeit within large relative errors. Between 0.8 fm and 1.4 fm, the lattice data lie noticeablybelow the pseudoscalar-octet exchange prediction.At distances up to 0.8 fm, one would certainly not expect the pseudoscalar-octet exchangeto saturate the integrand. We have attempted to model the higher-energy contributionsusing a constituent-quark loop, displayed in Fig. 7. For a constituent-quark mass of 350 MeV,the sum of this contribution and the pseudoscalar-octet exchange provides a good descriptionof the maximum height of the integrand. At distances | y | (cid:46) . . < | y | < . a type) would have the right sign to explainthe effect, since scalar-meson exchanges are known to contribute negatively to a hlbl µ . In thefuture, a calculation of the scalar contribution along the lines of [32] would be worthwhile tofind out whether it accounts for this missing effect. The charged-pion loop is also expectedto contribute negatively and we have studied this contribution in the scalar-QED frameworkin [18]. We found the integrand to be negligible compared with the π contribution beyond6 | y | = 0 . | y | around 1 fm.The effect of the finite volume on the lattice integrand is illustrated in Fig. 9. A cleareffect is seen between the comparatively ‘small’ ensemble H200 and the ‘large’ N202, theformer integrand lying below the latter. This finite-volume dependence matches in sign andtypical size the volume dependence of the pseudoscalar-octet exchange contribution, as seenin the figure. C. Lattice disconnected contribution
Fig. 12 (later in the text) displays the disconnected integrand for several gauge ensembles.It is compared to the prediction of the ( π , η ) exchange, including its appropriate weightfactor of − π , η ) exchange calculation is very small, and thelattice data from ensembles N202 and H200 confirms this expectation, at least at distances | y | < . π , η ) exchange already provides a decent description of the lattice datafor | y | below 1 fm; in other words, we do not observe a large short-distance contribution. The η (cid:48) and the σ mesons being singlets, they contribute to the disconnected diagrams with thesame weight factor as they would to the full HLbL amplitude [28, 30]. In order to estimatethe typical size of these heavy-meson contributions, we have computed the η (cid:48) contributionunder the following assumptions: The η (cid:48) mass was set to 982 MeV, a value close to the resultof our calculation on ensemble H101. Its coupling to two photons, F η (cid:48) γγ , was assumed tobe equal to its value at physical quark masses. The virtuality dependence of the transitionform factor can be modelled with a VMD ansatz, with vector mass 952 MeV. Under theseassumptions, the contribution is positive and sizeable (compared to the ( π , η ) exchange) upto | y | = 1 . σ meson,whose mass we have found to lie well below the ππ threshold. As a scalar, the σ wouldcontribute negatively and thus compensate to some extent the η (cid:48) contribution. Again, adedicated calculation in the framework of [32] would certainly be worthwhile. The estimatefor the η (cid:48) contribution displayed in Fig. 12 is only meant to be representative of one particularmeson-exchange contribution, with other cancellations being expected. V. RESULTS FOR THE QUARK-CONNECTED CONTRIBUTIONA. Results from Method 1
The lattice results for the quark-connected contribution to a hlbl µ using Method 1 have beengenerated along the direction y/a ∝ (1 , , , x is eitherlocal or point-split (conserved), while the three other currents are always local. Our resultsare summarized in Table V and the integrand for the ensemble N202 (also presented inSec. IV B, Fig. 7) is shown on Fig. 8. The signal-to-noise ratio clearly deteriorates rapidlyat large distances and we observe slightly better statistical precision when using a conserved7 − − . . . f ( | y | ) × [ f m − ] | y | [fm] localpoint-split . . . a c o nn µ × | y | max [fm] localpoint-split FIG. 8. Results for the ensemble N202 using Method 1 with Λ = 0 .
16. Black points use a localcurrent at the site x in Eq. (6) (which is our default approach) and blue points use a conservedcurrent at x . Left: integrand f ( | y | ) as a function of | y | . Right: integrated value of a conn µ as afunction of | y | max . vector current at x . Both discretizations give similar results, suggesting that there are smalldiscretization effects present. In the end we will quote a final result with the fully-localdiscretization for a direct comparison with the results obtained using Method 2.Although all the ensembles used in this work satisfy m π L >
4, we still expect significantfinite size effects (FSEs) due to the pseudoscalar-pole contribution, which is the dominantcontribution in model estimates of hadronic light-by-light scattering [33], being a long-range phenomenon. Even with heavy pions m π ≈
400 MeV, the tail extends beyond | y | =2 . | y | cut where the inte-grand is compatible with zero. For | y | > | y | cut , the tail is aproximated by the pseudoscalar-pole contribution in infinite volume. Finally, for | y | < | y | cut , the FSEs are estimated asthe difference between the pseudoscalar-pole’s contribution computed in finite and infinitevolume: a conn µ = a data µ + a tail µ + a FSE µ . (20)A systematic error of 25% is attributed to both the tail extension and the FSE correction.Since the same VMD model is used in both cases, we treat these corrections as being fullycorrelated.The value of a conn µ as a function of | y | cut is shown in the right panel of Fig. 9; we observe anice plateau for values | y | cut > . | y | cut . In particular, we note that the systematicerror on the FSEs is always larger than the statistical precision, except for our largestensemble N202.For the ensembles U103 and H200, with m π L <
5, we see very large FSE corrections, of8 − − . . . f ( | y | ) × [ f m − ] | y | [fm] H200N202 π + η ( L = 2 . fm) π + η ( L = 3 . fm) π + η ( L = ∞ ) . . . a c o nn µ × | y | cut [fm] raw lattice data (N202)FSE correctioncorrected lattice data FIG. 9. Study of FSEs for Method 1. Left: Integrand for the ensembles H200 and N202 with m π L = 4 . a conn µ for the ensemble N202 using the FSEcorrection prescription described in the text, as a function of | y | cut the order of 50 %. After these corrections the values of a conn µ do become compatible with theresults at m π L > σ . We also observe a systematic over-estimate in com-parison to the larger-volume results, and when it comes to our final continuum extrapolationwe will omit these results. O ( a )-improvement is not implemented for the vector currents used in this work, butour experience with other observables involving electromagnetic currents, such as the LOHVP [34] and the pion TFF [12], suggests the remaining O ( a ) terms are small comparedto the quadratic contribution. A linear fit in a leads to a conn , M1 µ = 104 . . × − with χ / d . o . f . = 0 .
4. To estimate the systematic error associated with this continuumextrapolation, we perform a constant fit which excludes the coarsest lattice spacing andobtain the slightly smaller value a conn , M1 µ = 96 . . × − with χ / d . o . f . = 0 .
3. Finally,we also tried a linear fit in the lattice spacing which leads to a conn , M1 µ = 113 . . × − .The results of all these fits are shown in Fig. 10.We quote our continuum-extrapolated value for the quark-connected contribution to a hlbl µ using Method 1 at the SU(3) f -symmetric point as a conn , M1 µ = 104 . . . × − , (21)where the first error includes both the statistical error and the systematic from the finite-size Label | y | cut [fm] a data µ a FSE µ a tail µ a conn µ U103 1.55 69.0(2.2) 31.0 4.3 104.3(2.2)(8.8)H101 1.73 71.6(2.9) 13.1 3.8 88.5(2.9)(4.2)B450 1.68 72.9(3.8) 18.3 4.0 95.2(3.8)(5.6)H200 1.54 71.1(2.8) 27.3 6.1 104.4(2.8)(8.4)N202 1.93 85.5(4.5) 9.0 1.7 96.3(4.5)(2.7)N300 1.59 79.2(4.2) 15.6 4.1 99.0(4.2)(4.9)TABLE V. Results for the connected contribution using Method 1 with four local vector currentsusing the decomposition given by Eq. (20). A 25% systematic to the total correction is used. . . . . a c o nn µ × a [fm ] fit linear in a fit linear in a constant fit FIG. 10. Continuum extrapolation of the connected contribution using Method 1. We perform alinear fit in a (red), a linear fit in a (green) or a constant fit (blue). The coarsest lattice spacingis excluded from the constant fit. Ensembles in grey have m π L < correction. The second error is an estimate of the continuum-limit extrapolation systematicerror, taken as half the difference between the linear in a and constant fit ans¨atze. B. Results from Method 2
In our measurement of the quark-connected contribution to a hlbl µ using Method 2 we focuson Λ = 0 .
4; this value was already indicated as being beneficial for the lepton loop as dis-cussed in Appendix B. We performed measurements on all ensembles with Λ = 0 . , . , . , and 1 . . . . x and/or z . We found that putting a conserved current at z yields a result roughlycomparable (point-by-point) to the determination with just 4 local currents. Having aconserved current at x appears to introduce large, unwanted discretization effects. Froma computational standpoint the calculation with four local currents is simpler and has noapparent downside, so that is what we will present from here on.The left plot of Fig. 11 illustrates the finite-size effect between ensembles N202 and H200,and the discrepancy between these ensembles that differ only by volume is significant. Thepion-pole prediction describes the tails of both data-sets reasonably well at large-enoughvalues of | y | , although it completely under-estimates the position and height of the short-distance peak of the integrand. It is also worth noting how less statistically precise the resultof H200 is compared to N202 for a comparable number of measurements; this gives someindication that the statistical precision is linked to either m π L or the physical volume. Itis clear that much like for Method 1 there is a significant signal-to-noise problem for largevalues of | y | .0 Label | y | cut [fm] a data µ a FSE µ a tail µ a conn µ U103 1.85 58.9(1.8) 14.3 13.3 86.5(1.8)(6.9)H101 2.55 75.7(2.9) 6.9 2.6 85.2(2.9)(2.4)B450 2.00 73.6(3.3) 7.6 9.4 90.3(3.3)(4.3)H200 1.75 68.6(1.8) 13.7 13.6 95.8(1.8)(6.8)N202 2.60 91.0(2.5) 3.8 1.9 96.7(2.5)(1.5)N300 2.10 79.0(1.8) 7.1 6.2 92.3(1.8)(3.4)TABLE VI. Results for the connected contribution using Method 2. Here a data µ corresponds to thevalue using lattice data up to some linearly-interpolated value of y = | y | cut , chosen to minimize thetotal error of a conn µ . Again, a 25% systematic to the total correction is used. In the last column wegive the infinite-volume result. − − f ( | y | ) × [ f m − ] | y | [fm] H200N202 π + η ( L = 2 . fm) π + η ( L = 3 . fm) π + η ( L = ∞ ) a c o nn µ × | y | cut [fm] raw lattice data (N202)FSE correctioncorrected lattice data FIG. 11. Study of FSEs for Method 2. Left: Integrand for the ensembles H200 and N202 with m π L = 4 . . a conn µ for the ensemble N202 using the FSEcorrection prescription described in the text, as a function of | y | cut . We perform the same finite-size correction procedure for Method 2 as we did for Method1 above. On the right of Fig. 11 we show the stability of performing the FSE correction withvaried | y | cut matching point on ensemble N202. We find excellent stability for many differentvalues of | y | cut and the matching point of 2 . a conn,M2 µ = 98 . . × − . (22)1 Calculation N Conf N Src N Prop
Method 1 800 48 280 , , , N Conf gives the number ofgauge configurations used, N Src the number of source positions per gauge configuration, and N Prop the total number of propagator solves.
The quoted error is a combination of the statistical and 100%-correlated finite-size system-atic. We observe that within error this value is in complete agreement with the determinationof Method 1.
C. Comparison of the two methods
The appeal of Method 2 for computing the connected contribution is mostly practical: itis computationally far less expensive than Method 1. The saving between the two is roughlyan order of magnitude, see Tab. VII for a exemplary comparison of computational cost forone particular ensemble, N202. This is because Method 2 effectively replaces sequentialpropagator solves by additional, much cheaper, QED kernel evaluations. The downsideof using these additional kernels is that their combination tends to broaden the integrand f ( | y | ). This behavior was seen in the lepton loop study of Appendix B and is also clearlythe case with all the lattice QCD data in this work. We can use the parameter Λ of Eq. (7)to partially ameliorate this broadening; the use of such a subtraction kernel appears to bevery important specifically for Method 2, as this regulator offers little to no advantage forMethod 1.If we compare the results for the integrand of H200 and N202 using Method 1 with thoseof Method 2 (Figs. 8 and 11 respectively) we see that the integrand for Method 2 is ingeneral less-peaked at short distances and extends further in | y | . For example, the integrandfor N202 using Method 1 is effectively zero around 2 fm, whereas for Method 2 it becomeszero closer to 3 fm. This behavior is reflected in Tables V and VI by larger choices of | y | cut for Method 2 compared to Method 1.Again comparing Tables V and VI, we see that for both methods the smaller boxes (thosewith m π L <
5) require a significant finite volume correction. As the data for Method 2 usesthe direction (1 , , ,
3) we approach the boundary of the lattice ( L/
2) slowly for increasing | y | , and so the finite volume effect is smaller in comparison to the direction used in Method 1.We do, however, see somewhat larger discretization effects for Method 2 compared to thosefound in Method 1, perhaps O (15%) opposed to O (10%) at our coarsest lattice spacingrespectively. It is quite possible this is due to the Λ-regulator enhancing the integrand atshorter distances. Nevertheless, this is not a significant problem as we have several finelattice spacings to help determine the continuum limit.If we were to use Λ = 0 . m π L , but this wouldbecome much more problematic for lighter-pion-mass ensembles where the signal is expectedto degrade quickly at large distances and the integrand is expected to be even broader.2In the following sections, when we combine the results for the connected and discon-nected contributions, we will use the results from Method 2 as our connected contribution.This is because they are statistically more precise while still being consistent with those ofMethod 1. VI. RESULTS FOR THE QUARK-DISCONNECTED CONTRIBUTION -100-50 0 50 100 150 0 0.5 1 1.5 2 2.5 f( | y | ) x [f m - ] |y| [fm] kernel L ( Λ ) , Λ =0.4(m π , m V , f π )/MeV = (416, 924, 105.93)(m η ’ , m VMD , f η ’ )/MeV = (982, 952, 73.6)L H101 =2.76fm, L
U103 =2.07fm
H101 U103 dir=(1,1,1,3)Conn: π o + η inf.volDisc: π o + η inf.volDisc: π o + η L=L
H101 dir=(2,2,2,0)Disc: π o + η L=L
H101 dir=(1,1,1,3)Disc: π o + η + η ’ inf.vol -100-50 0 50 100 150 0 0.5 1 1.5 2 2.5 f( | y | ) x [f m - ] |y| [fm] L N202 =3.08fm, L
H200 =2.06fm(m π , m V , f π )/MeV = (410.56, 952, 111.59)(m η ’ , m VMD , f η ’ )/MeV = (982, 952, 73.6)kernel L ( Λ ) , Λ =0.4 N202H200 dir=(1,1,1,3)Conn: π + η inf.volDisc: π + η inf.volDisc: π + η L=L
H200 dir=(2,2,2,0)Disc: π + η L=L
H200 dir=(1,1,1,3)Disc: π + η + η ’ FIG. 12. Integrand for the disconnected contribution on ensembles H101 and U103 (left), as wellas N202 and H200 (right). The lattice data are shown as black points. The black dashed line showsthe fully-connected model prediction, the blue line the π + η contribution for the disconnectedand the magenta line gives the π + η + η (cid:48) . The blue points show the π + η contribution in finitevolume. Table VIII lists the ensembles and statistics used for the computation of the quark-disconnected contribution. As the smaller ensembles (U103, H101, B450, H200) were con-siderably cheaper to perform inversions on, their statistics is greatly enhanced. As the lat-tice volume increases, the cost of propagator inversions increases with some power V n , with n >
1, and this quickly becomes the dominant cost of the computation. The column N Src indicates the number of point-source propagators inverted per gauge configuration to buildthe grid and the final column indicates the maximum and minimum number of equivalentvalues of | y | available from the set-up. For the ensembles with open boundary conditions,the number of self-averages, N Equiv , for a given | y | decreases as | y | /a increases. Thereforeopen temporal boundaries make this calculation much more difficult as the signal degradesrapidly with large | y | /a .The ensemble B450 has lattice volume 32 ×
64 and periodic boundary condition in time,so a fully-periodic grid built of the two basis vectors (4 , , ,
0) and (2 , , ,
4) was used. Allof the other determinations had multiples of (1 , , ,
3) and (2 , , , | y | directions, so it ispossible that B450 might have noticeably different discretization and finite-volume effects asthis direction lies in a different lattice irreducible representation of the broken rotation group O (4). However, as we see in Figs. 12 and 13 the short-distance contribution to the integral issmall, and so we can assume the same is true of the discretization effects. As we appropriatelycorrect for FSEs with this choice of direction, we do not expect a significant discrepancycompared to the open-boundary data. Upon continuum extrapolation (Fig. 14) it does seem3 Label ( t min , t max ) /a N Src N Conf N Equiv (max,min)U103 (20 , ,
72) 272 1008 (1568,128)B450 ∗ (0 ,
64) 128 1611 (512,128)H200 (24 ,
72) 272 500 (1568,128)N202 (36 ,
93) 480 225 (2784,672)N300 (27 ,
99) 600 271 (3504,1152)TABLE VIII. The setup used for each ensemble in the computation of the quark-disconnectedcontribution. N Src gives the number of propagator inversions per configuration, N Conf gives thenumber of gauge configurations used and N Equiv gives the maximum and minimum number ofequivalent values of | y | averaged per configuration. Shorter separations have larger numbers ofself-averages, whereas larger values of | y | have smaller N Equiv . The ensemble B450 has a periodictemporal boundary and temporal length of 2 × that of the the spatial, so values of y/a = n (2 , , , that this ensemble is consistent with the others, indicating that rotation-breaking artifactsare not the main source of discretization effects.The integrand for the disconnected contribution is displayed in Fig. 12 for the valueΛ = 0 .
4; much like for Method 2 we find this value to be preferable. The integrand iscompared to the prediction for the exchange of the π and η mesons with a VMD TFF.In addition, the same prediction including an estimate of the η (cid:48) contribution, based on theassumptions in Section IV C, is indicated.We do not see any statistically significant finite-volume effects in the integrands betweenthe ensembles U103 and H101, and H200 and N202, for | y | < π , η ) exchange in finite volume. The central valuesof the integrand obtained on ensembles with different volumes differ substantially at somelarger values of | y | , however this is also in the regime where the signal is rapidly deteriorating,if not lost already. There is a trend in the tail for the larger-volume results to enhancethe magnitude of the disconnected contribution, much like what we saw in the connectedcontribution. We see some enhancement of the integrand compared to the π + η + η (cid:48) prediction at short distances; the likely cause of this is the contribution from scalar mesons.Clearly, the π and η exchanges already provide a rather good description of the shapeof the integrand, unlike in the case of the connected contribution (e.g. Figs. 8 and 11). Atthe same time, the loss of the signal beyond 1.2 fm means that we cannot, at this point,confirm the validity of the ( π , η ) exchange at long distances. Given the rapid degradationin signal of the disconnected data, probing large distances of the disconnected contributionwill be a very challenging undertaking. There are however good reasons to believe that thisdescription should apply in that regime, and in the following we assume this to be the case.Table IX summarizes our results. We perform the FSE matching at a single linearly-interpolated point of | y | cut = 1 . π , η ) exchange.This correction is of the order of 100% of the lattice-determined contribution. We have also4 − − − − − . . . a d i s c µ × | y | cut [fm] raw lattice dataFSE correctioncorrected lattice data − − − − − . . . a d i s c µ × | y | cut [fm] raw lattice dataFSE correctioncorrected lattice data FIG. 13. Value of a disc µ as a function of | y | cut for ensembles H101 (left) and N202 (right).Label a data µ a FSE µ a tail; π + ηµ a tail; η (cid:48) µ a disc µ U103 − . . − . − . − . . . − . . − . − . − . . . − . . − . − . − . . . − . . − . − . − . . . − . . − . − . − . . . − . . − . − . − . . . a hlbl µ at the SU(3) f -symmetricpoint. A breakdown of the FSE contributions to the total results are shown, again a 25% systematicto the total finite-size correction is used. The value of | y | cut was 1 . computed an estimate of the η (cid:48) exchange contribution as described in section IV C. Thevalues in Table IX show that its contribution to the tail, | y | > . η (cid:48) to the tail in the centralvalue of a disc µ .Anticipating the combined analysis presented in section VII A, the continuum-extrapolateddisconnected value we obtain from a constrained-slope fit to both the connected Method 2data and the disconnected data is a disc µ = − . . × − . (23)We observe that this result amounts to be about ( − /
3) of the connected contribution.
VII. THE TOTAL a hlbl µ The main purpose of this section is to describe how we arrive at our final result for thetotal a hlbl µ (subsection VII A), in which the systematics of the continuum extrapolation arediscussed in detail. The following sub-section VII B presents a study of the dependence ofour result on the muon mass, and finally, subsection VII C discusses what outcome one mayexpect for a hlbl µ at physical quark masses based on our findings.5 A. Final result for a hlbl µ in SU(3) f -symmetric QCD
1. Combining the connected and disconnected contributions
In order to obtain the final result for a hlbl µ = a conn µ + a disc µ , we combine the disconnecteddata with the connected data of Method 2, as it is consistent with, yet statistically moreprecise than that of Method 1. We then perform two analyses on this combined data set.The full set of data for the two contributions and their sum is shown in Fig. 14, along withfits and results from the first analysis whose description follows. Both the connected anddisconnected data show a negative slope in a despite the two contributions having oppositesigns. This indicates that the leading discretization effects do not arise from a commonmultiplicative effect such as the renormalization factor.The first analysis consists in simultaneously extrapolating both the connected and thedisconnected contributions to the continuum with a slope in a constrained to be equal forthe two contributions, before summing the two constant parameters of the fit to obtain a hlbl µ .In this procedure, the results for the connected contribution and the disconnected contri-bution are both corrected for finite-size effects and for the extension of the | y | -integrand.The statistical errors are obtained under the bootstrap, and a correlated gaussian samplewith appropriate width is propagated for the systematic error due to the correction. Thesystematic error of the correction, which is set to 25% of its size, is treated as being fullyanti-correlated between the connected and the disconnected contribution (recall that thetwo contributions have opposite sign). The fit ansatz for the two continuum extrapolationsis a polynomial of degree one in a . For this analysis we obtain the result, a hlbl µ = 65 . . × − . (24)In the second analysis method, the connected and disconnected contributions are summedprior to performing the continuum extrapolation. This extrapolation is then performed lin-early in a . In this analysis, the correction applied to the lattice data is split into two parts,(a) a FSE µ , the pure finite-size correction on the integrand, and (b) a tail µ , the extension of theintegrand beyond some | y | cut . In combining the connected and disconnected contributions,each part is treated as being fully anti-correlated between the connected and the discon-nected contribution; however, in contrast with the first analysis, the correlation between thesystematic errors of the two parts is considered to be zero. The alternative point of viewadopted here is that predicting the tail of the infinite-volume integrand is a rather differentexercise from predicting the finite-size effect on the finite-volume integrand; one of the twopredictions could be successful while the other is not. The result of this procedure is a hlbl µ = 64 . . × − . (25)The two analyses produce compatible values for a hlbl µ but the first is more precise andallows more flexibility since the two contributions are treated separately. Therefore, we willuse the first analysis for the final result and study variations on it to estimate the systematicuncertainty.
2. Continuum extrapolation systematics
Our data are noisy and it is difficult to find a fit that describes them perfectly, so weidentify several different choices of continuum extrapolation to investigate the spread and6 a [fm ] -60-40-20020406080100 a µ h l b l x ConnectedConn. + Disc.Disconnected
FIG. 14. Combined continuum-extrapolation analysis for the connected and disconnected data.The fits were constrained to have the same slope and the final sum of the two was performed onthe constant fit parameters. Also shown is the individual sum of the connected and disconnectedpieces. The purple cross represents the addition of the continuum-extrapolated results for theconnected and disconnected contributions. provide an associated continuum-extrapolation systematic. We note that (also expressed inSec. V A) as we are not using the O ( a )-improved vector currents it is possible that we havea term linear in the lattice spacing. We consider the two fit forms: a conn µ ( a ) = a conn µ (0) + Aa n , a disc µ ( a ) = a disc µ (0) + Ba m , (26)with various cuts to the data, constraints on A and B , and choices of n and m as listed inTable X. We then add the distributions a hlbl µ = a conn µ (0) + a disc µ (0) to obtain the results shownin Fig. 15. Index a -Cut [fm] Constraint n m Index a -Cut [fm] Constraint n m A = B A = B B = 0 1 - 8 None B = 0 2 -4 < . < . < . A = B < . A = B It is clear that any time we omit the coarse data the fit wants to flatten the slope. This isdue to the anomalously low result of N300. When we do not include the coarse ensembles theerror increases substantially, which is fairly indicative that the fit is struggling to accuratelymodel the data.7
Fit Index a µ h l b l x Estimated SystematicFit Result
FIG. 15. An estimate of our continuum extrapolation systematic on the full a hlbl µ for the combina-tion of Method 2 and the disconnected data. If we perform a fit linear in a the central value moves up, as was also the case forthe Method 1 continuum extrapolation. We choose to quote the linear in a fit to all thedata with a constrained slope as our final result as it has the best χ / d . o . f . = 2. For thecontinuum-extrapolation systematic, we use the lower error bound of the largest fit resultand the upper bound of our lowest fit result. It is clear that constraining the slope or lettingit vary does little to the position of the central value apart from reducing the error, thissuggests that the fit is having a hard time accurately determining the slope with the qualityof the data we have at present. B. The dependence of our results on the muon mass
Since the mesons at the SU(3) f -symmetric point can be viewed as ‘heavy’ degrees offreedom relative to the muon, we would expect a hlbl µ to be roughly proportional to m µ . Here,we will study what happens if we re-scale the muon mass on one of our ensembles. We aremotivated to do so by our experience on related projects [34] whereby adjusting the muonmass by some dimensionless ratio can flatten the approach to the chiral limit. Here weinvestigate this idea by defining two quantities, a hlbl µ = (cid:18) f Latt. π f Phys. π (cid:19) a hlbl µ ( m Phys. µ ) , (27)and, ˜ a hlbl µ = a hlbl µ (cid:18) f Latt. π f Phys. π m Phys. µ (cid:19) . (28)The first quantity rescales the integrated result by the lattice-determined pion decay constantdivided by the value in continuum, squared. The second quantity re-scales the muon mass8 | y | max [fm] a µ c onn ( | y | m a x ) a µ conn ( f π Lat /f π Phys. ) m µ ( f π Lat /f π Phys. )Not rescaled
FIG. 16. Partially-integrated lattice results for the ensemble H101 with and without re-scaling ofthe integral. used as input in our determination by this ratio. These two definitions would be comparableif a hlbl µ scales as m µ , which is to be expected in the heavy-quark limit.Fig. 16 illustrates the effect of these re-scaling procedures on the connected contributionto a hlbl µ on one of our coarsest and largest boxes, H101 (results are from Method 2 withΛ = 0 . very naively infer how much we expect the result to grow as weapproach the chiral limit, and it appears that for the connected contribution this could beof the order of a 25% increase. C. Expectations for a hlbl µ in QCD with physical quark masses In this subsection, we first quantify the contributions to a hlbl µ at the SU(3) f -symmetricpoint not coming from the light pseudoscalars. These contributions are expected to have onlya modest quark-mass dependence, and they are also the hardest to determine quantitativelyby (experimental-)data-driven methods. Therefore any information on these contributionsis worthwhile collecting. By using our determination of this contribution together with ourprevious calculation of the π exchange, we can arrive at an estimate for a hlbl µ at physicalquark masses.We begin by noting that, subtracting the π and η contributions (respectively Eqs. (16)and (18)) from our final SU(3)-point result (31), the contribution of heavier intermediatestates amounts to a hlbl , SU(3) f µ − a hlbl ,π + η, SU(3) f µ = (37 . ± . × − . (29)In particular, this contribution accounts for 57% of the total a hlbl µ , and we have added all9statistical and systematic errors in quadrature.Next, we may try to roughly estimate a hlbl µ at the physical point. As the ( u, d ) quarkmasses are lowered to their physical values at fixed trace of the quark mass matrix, it is thepion whose mass changes by the largest factor: it becomes a factor of three lighter. Sincewe have an evaluation of the π exchange contribution at the SU(3) f -symmetric point andat the physical point (see Eqs. (16–17)), we can correct for this effect, a hlbl , SU(3) f µ − a hlbl ,π , SU(3) f µ + a hlbl ,π , phys µ = (104 . ± . × − . (30)One can think of Eq. (30) as a rough estimate of a hlbl µ at the physical point, purely basedon our lattice QCD results and the assumption of a negligible quark-mass dependence ofthe non- π -exchange contributions. As argued in the next paragraph, we expect such anestimate to hold at the 20% level. It is well in line with the most recent evaluations [5].In order to assess the systematic uncertainty of such an estimate, we perform a slightlymore sophisticated method to correct for the quark-mass dependence of the η contributionand the charged pion loop. Within the scalar QED framework, we find − . × − forthe pion loop, to be doubled to include the kaon loop, and we expect a factor of two tothree reduction if one includes an electromagnetic form factor for the pseudoscalar. There-fore, further subtracting a hlbl , ( π ± ,K ± ) , SU(3) f µ ≈ × ( − . × − from Eq. (29), we obtain(43 . ± . × − . This number represents our estimate of the non-pseudo-Goldstonecontributions at the SU(3) f -symmetric point . We observe that by neglecting the quark-mass dependence of this contribution, and using the dispersive π -exchange [5, 35] and theCanberbury-approximant η -exchange [5, 26] results for a hlbl ,π + η, phys µ = 79 . × − andthe box contribution [5, 36] to a hlbl , ( π ± ,K ± ) , phys µ = − . × − , we arrive at the estimate10 a hlbl , phys µ = 79 . . − . . .
3) = 106 . . a hlbl µ at the physical point. This uncertainty estimate alsogenerously covers the a hlbl , phys µ value obtained by assuming that the non-pseudo-Goldstonecontributions increase from the SU(3) f to the physical point by a factor ( f SU(3) f π /f phys π ) to account approximately for the quark-mass dependence of the QCD resonances; see theprevious subsection concerning this point. VIII. CONCLUSIONS
In this work we have computed the hadronic light-by-light contribution to the g − f -symmetric point with m π = m K ≈
420 MeV. Wechose to initially work at the symmetric point for several reasons: Due to the significantlyreduced computational cost (as compared to simulations at physical quark masses), weare able to control all known sources of systematic error, in particular the finite-size andfinite lattice spacing effects. Second, the SU(3) f -symmetry implies that only two out of fivequark-contraction topologies contribute, and it simplifies the hadronic models with whichthe integrand can be compared and interpreted. For instance, the η -exchange contributionsimply amounts to one third of the pion-exchange. Third, the overall contribution from Note that since the lightest f meson is stable, there is no reason to treat its contribution as part of the ππ rescattering effect. a hlbl µ , we have performed an exploratory studyof the low-lying meson spectrum at the SU(3) f point. The most remarkable feature isthe existence of a stable singlet J P C = 0 ++ meson with a mass of about 680 MeV. Also,our previous calculation of the pion transition form factor [12] allows us to quantify thecontributions of the π and η exchanges.Our strategy to calculate a hlbl µ relies on coordinate-space perturbation theory, for whichmuon and photon propagators are computed in infinite volume. We have presented theintegrand for the final, one-dimensional integral over the distance | y | of a quark-photon ver-tex from one of the other two internal vertices, since it contains more information than thefinal a hlbl µ value. For the quark-connected contribution we have identified two methods ofcalculation that we call Method 1 and Method 2. The former amounts to a direct computa-tion of the three connected diagrams, but it is a numerically costly approach, as it involvesmany sequential-propagator calculations. To make this computation much cheaper, we haveutilized several changes of variables and translational invariance to rewrite the integral interms of an easy-to-calculate single diagram and a combination of different kernels; this wecall Method 2.For a single | y | -value, Method 2 requires only two propagator inversions, at the meagrecost of a more-complicated QED-kernel calculation. Method 2 has another computationaladvantage over Method 1 in that it allows one to store a set of propagators in memory andperform their integrals, redefining the origin to be each propagator source; thereby, for N propagators we can compute N ( N −
1) non-zero samples of f ( | y | ). If the source points areevenly spaced and the volume is periodic this amounts to N self-averages per | y | . Suchself-averaging is crucial in reducing the cost of the calculation.We note that the combination of kernels needed for Method 2 broadens the integrandin | y | significantly. To counteract this effect, we have introduced a new class of subtractedkernels with a gaussian-regulator Λ. This parameter Λ effectively allows us to tune theshape of the integrand to peak at shorter-distances. We find that a value of Λ = 0 . π and η meson exchange contributionto the integrand using a vector-dominance transition form factor and compared it to thelattice data; only at fairly large | y | does the prediction quantitatively represent the inte-grand. We therefore attempt to conservatively incorporate as much lattice data as possiblebefore making contact with the model. For the disconnected contribution, the model doesa satisfactory job of describing the data over the entire range of | y | where we have signal,and we therefore match on to the prediction at shorter distances, where we still have controlover the statistical errors of the lattice data.1We have shown that both the disconnected and connected contributions have a non-negligible discretization effect within our measured precision. It appears to be of the samesign and comparable magnitude for both contributions, but ultimately this extrapolationappears well under control. We decide to quote a result from a fit linear in a with aconstrained slope to the data of Method 2 and the disconnected as, a hlbl µ = (65 . ± . ± . × − , (31)at the SU(3) flavor-symmetric point, where the first error results from the uncertainties onthe individual gauge ensembles, and the second is the systematic error of the continuumextrapolation.We have discussed in subsection VII C how we expect our result for a hlbl µ to evolve as theup and down quark masses are lowered towards their physical values at fixed trace of thequark mass matrix. Correcting for the increase in the π exchange contribution using ourprevious lattice calculation [12], we arrive at a value (Eq. (30)) which is very well in linewith the most recent phenomenological [5] and lattice QCD results [16]. This value is quitestable under varying the assumptions about the quark-mass dependence of heavier-statecontributions.In order to reduce the systematic uncertainty of a hlbl µ at physical quark masses usinglattice QCD, obviously simulations at lighter quark masses are needed. The methods wehave developed to correct for finite-size effects and to extend the tail of the | y | -integrandbased on the π exchange will be extremely valuable in this endeavor. While we have reacheda semi-quantitative understanding of the integrand in terms of hadronic models, further workis needed on the theory side to bring this description to a fully quantitative level. ACKNOWLEDGMENTS
This work is supported by the European Research Council (ERC) under the EuropeanUnion’s Horizon 2020 research and innovation programme through grant agreement 771971-SIMDAMA, as well as by the Deutsche Forschungsgemeinschaft (DFG) through the Col-laborative Research Centre 1044 and through the Cluster of Excellence
Precision Physics,Fundamental Interactions, and Structure of Matter (PRISMA+ EXC 2118/1) within theGerman Excellence Strategy (Project ID 39083149). The project leading to this publica-tion has also received funding from the Excellence Initiative of Aix-Marseille University -A*MIDEX, a French “Investissements d’Avenir” programme, AMX-18-ACE-005. Calcula-tions for this project were partly performed on the HPC clusters “Clover” and “HIMster II”at the Helmholtz-Institut Mainz and “Mogon II” at JGU Mainz. Our programs use the de-flated SAP+GCR solver from the openQCD package [37], as well as the QDP++ library [38].The spectroscopy calculation also made use of the PRIMME library [39] and was partly per-formed on the supercomputer JUQUEEN [40] at J¨ulich Supercomputing Centre (JSC); weacknowledge the Gauss Centre for Supercomputing e.V. ( )for providing computer time. We are grateful to our colleagues in the CLS initiative forsharing ensembles. For this purpose, a model-independent parametrization of the pion transition form factor is used, asopposed to the simple vector-meson dominance parametrization. Appendix A: Pseudoscalar-meson-exchange-channel matching in Method 2 and thequark-disconnected contribution
In this work, when we compute the quark-connected diagrams using Method 2 and thequark-disconnected diagrams, only a subset of Feynman diagrams are considered due to thecontributing Wick-contractions at QCD level. In order to study the FSE in an effectivefield theory framework, one has to map these contractions to the corresponding Feynmandiagrams in the effective theory. Our study of the FSE is based on pseudoscalar-meson(PS) exchange from Chiral Perturbation Theory. The matching between the QCD Feynmandiagrams and the PS-meson-exchange channels can be done in various ways. Here, we presenta determination of the matching using Partially-Quenched QCD (PQQCD) and Partially-Quenched Chiral Perturbation Theory (PQChPT). There have been many applications ofPQQCD in the Lattice QCD community (see e.g. [41–43]). In particular, it can serve asa tool to give an estimate of the size of the contribution coming from quark-disconnecteddiagrams compared to the connected ones for some observables [44]. Here we do not intendto detail the formalism, but only to give some arguments to explain how we reach themappings explained in our methodology.PQQCD is a theory with a graded Lie-group SU( N | M ) as symmetry group. In thistheory there are N − M sea quarks and M valence (quenched) quarks with their ghostcounterparts. The presence of these quenched quarks does not change the partition functionfrom the un-quenched theory with N − M dynamical quarks. Therefore, one can formulateeach of the different Wick contractions needed for computing an n -point function (in anun-quenched theory) in a partially-quenched theory by adding certain number of quenchedquarks. Then, one can study the long-distance behavior of the observable using PQChPT,which is the corresponding effective field theory to PQQCD.With only three flavors (as is the case in this work), one cannot build a four-pointfunction that requires only quark-connected diagrams, because each quark line would requirea different flavor in order to avoid disconnected diagrams. However, such a four-pointfunction can be constructed with an additional quenched quark flavor, i.e. SU(4 |
1) PQQCD.For instance, if we introduce a quenched quark r and its ghost ˜ r with the same mass as theones originally present in the SU(3) f theory, namely u , d and s , in the computation of thefour-point function (cid:104) ¯ uγ µ d ( x ) ¯ dγ ν s ( y )¯ sγ σ r ( z )¯ rγ σ u (0) + h.c. (cid:105) , (A1)only quark-connected diagrams are involved. Eq. (A1) is exactly the four-point function (cid:101) Π (1) µνσλ ( x, y, z,
0) that is used for the connected contribution in Method 2.Moreover, one can study the quark-disconnected diagrams that we consider in this workwithin the same partially quenched theory. For instance, (cid:104) ¯ uγ µ d ( x ) ¯ dγ ν u ( y )¯ sγ σ r ( z )¯ rγ σ s (0) + h.c. (cid:105) (A2)only contains the disconnected diagram where ( x, y ) and ( z,
0) are connected by quark lines.To determine the relevent PS-meson-exchange channel, one can consider the Wess-Zumino-Witten term [45, 46] in PQChPT at leading order in the decay constant. Eq. (A1)can be understood as a four-point vector current correlation function. A general procedureis to include all the possible super traces of operators [47]. However, in the case of the cou-pling to external charged vector currents, which are the Noether currents of the symmetrygroup SU( N | M ) whose generators are super-traceless, the only allowed non-vanishing term3 + ∼∼ − FIG. 17. Correspondence between connected (Method 2) and disconnected diagrams (left) and PS-meson exchange channels (right) in position space of the four-point functions described in the text.The top-left figure should be understood as containing quark flow in both orientations, representingexactly all the relevent diagrams for Method 2. Note that the charge factors are not included here. is proportional to str( T a T b T c ) (cid:90) d xφ c (cid:15) µνρσ F aµν F bρσ , (A3)where φ is the Goldstone boson/ghost field, F µν is the field strength of the external vectorcurrent and the T a ’s are the generators of the group. As a result, in momentum-space,the PS-meson exchange in the s -channel of the two-vector-currents-to-two-vector-currentsprocess v a v b → v a v b is proportional tostr( { T a , T b } T c )str( { T a , T b } T c ) g c c G ( p ) , (A4)where G is the scalar field propagator , G ( p ) = 1 p + m π , (A5)and g ab = 2str( T a T b ) . (A6)The matching between different contractions and the PS-exchange channels based on thiscomputation is shown in Fig. 17. With the charge factors included, one will recover the ratiobetween the total quark-disconnected contribution and the total quark-connected given inRef. [28], based on charge factor arguments. Appendix B: The new subtraction with the lepton-loop
In this section we investigate the new subtraction for both Methods 1 and 2 using theinfinite volume lepton-loop. A comparison of previous results with our kernel to the lepton-loop can be found in [17, 27]. The lepton-loop result can be determined easily by standardnumerical integration techniques, the results of which can be used to inform our choices ofkernel parameters in the full QCD calculation.For both figures (18 and 19) we set the lepton mass equal to that of the muon.Fig. 18 shows the results for Method 1. There is little difference between the choicesof Λ for this calculation and they all saturate the integral reasonably quickly. The result4when Λ = ∞ is too peaked at short distances to make this a useful kernel for the latticecomputation. All of the peak positions for the integrands are roughly at the same (small)value of | y | .In Fig. 19 we show the results for Method 2. We see that the usual kernel (Λ = ∞ ) is verypeaked at short distances, this will make it unsuitable for the lattice calculation. We alsonote that Λ = 0 gives a broad integrand that very slowly approaches the exact result, thisalso makes it somewhat unsuitable for a lattice calculation unless a large physical volumeand large values of | y | are available.The choice of Λ = 0 . | y | than Λ = 0and it goes to zero at large distances rapidly. There is a small negative contribution to theintegral at large | y | , but in the full QCD case can be corrected for by modelling the tail, as Λincreases this negative correction grows. If our lepton-loop result translates to full QCD thischoice of Λ should allow for most of our integral to be encompassed by the lattice volume. |y| / m µ -1e-0901e-092e-093e-094e-095e-096e-097e-098e-099e-091e-081.1e-081.2e-08 f( | y | ) m µ Λ = 0 Λ = 0.4 Λ = 0.8 Λ = 1.0 Λ = ∞ (a) Method 1 integrands |y| max / m µ a µ L e p . ( | y | m a x ) Λ = 0.0 Λ = 0.4 Λ = 0.8 Λ = 1.0 Λ = ∞ Exact Result (b) Method 1 integrals
FIG. 18. Left: the integrands in Method 1 of several Λ choices of the new subtracted kernel.Right: the results of the partial integration of the integrands up to some | y | max . These results wereobtained with lepton mass equal to the muon mass and the exact result used was 4 . × − [48]. If we compare the two methods using the lepton-loop we see that the parameter Λ makesa significant difference to the integrand for Method 2. Any of the choices Λ = 0 . , . , . , . |y| / m µ -1e-0901e-092e-093e-094e-095e-096e-097e-098e-099e-091e-081.1e-081.2e-08 f( | y | ) m µ Λ = 0.0 Λ = 0.4 Λ = 0.8 Λ = 1.0 Λ = ∞ (a) Method 2 integrands |y| max / m µ a µ L e p . ( | y | m a x ) Λ = 0.0 Λ = 0.4 Λ = 0.8 Λ = 1.0 Λ = ∞ Exact Result (b) Method 2 integrals
FIG. 19. Same as Fig. 18 but for Method 2. seem usable for Method 1. It is interesting to note that Method 2 actively shifts the peak ofthe integrand to typically intermediate distances. So, with a careful choice of the parameterΛ we can tune the kernel to peak in a region where the calculation is not so sensitive toeither of the lattice cut-offs of lattice spacing and volume.
Appendix C: Kernel independence of a conn µ and a disc µ The validity of the subtractions applied to the QED kernel originates from the Wardidentity (current conservation) satisfied by the electromagnetic current. Thus, one willexpect a kernel-independent result at the level of the total result of a hlbl µ . In this paper,we have treated the contribution from the connected diagrams ( a conn µ ) and the contributionfrom the disconnected diagrams ( a disc) µ ) separately. One natural question would be whetherthese two quantities themselves are also kernel-independent. The answer to this question isyes, due to the Ward identity satisfied by the vector current. We will give a proof of thisstatement in the SU(3) f case.In Euclidean space-time, if an infinitesimal transformation δ (cid:15) of the fields does not change6the path integral measure, we have (cid:104) δ (cid:15) S O(cid:105) = (cid:104) δ (cid:15) O(cid:105) , (C1)where S is the QCD action and O is any operator. Without changing the action, we canconsider a partially-quenched theory by introducing a quenched quark r and its ghost partner˜ r (see Appendix A). If we consider an infinitesimal transformation generated by the localSU(2) f symmetry of the quark pair ( u, d ) and their conjugates : δ (cid:15) u ( x ) = i(cid:15) ( x ) d ( x ) , δ (cid:15) d ( x ) = i(cid:15) ( x ) u ( x ) ,δ (cid:15) ¯ u ( x ) = − i(cid:15) ( x ) ¯ d ( x ) , δ (cid:15) ¯ d ( x ) = − i(cid:15) ( x )¯ u ( x ) , (C2)we have the following expression for the variation of the action related to the correspondingNoether current δ (cid:15) S = i(cid:15) ( x ) ∂ λ (¯ uγ λ d + ¯ dγ λ u )( x ) . (C3)We now consider the trilocal operator O µνσ ( x, y, z ) := ¯ uγ µ d ( x )¯ sγ ν r ( y )¯ rγ σ s ( z ), which trans-forms as δ (cid:15) O µνσ ( x ) = i(cid:15) (¯ uγ µ u − ¯ dγ µ d )( x )¯ sγ ν r ( y )¯ rγ λ s ( z ) . (C4)Plugging in everything in Eq. (C1), we have ∂ λ (cid:104) ¯ uγ µ d ( x ) ¯ dγ λ u ( x )¯ sγ ν r ( y )¯ rγ σ s ( z ) (cid:105) = ∂ λ (cid:104) Π µλ ( x, x )Π νσ ( y, z ) (cid:105) U = 0 , (C5)where Π µν ( x, x )Π νσ ( y, z ) gives one of the contractions needed for the disconnected compu-tation. One can then easily show that all the disconnected diagrams have similar current-conservation property. Consequently, the subtractions that we apply to the QED-kernel donot change the result on the contribution from the disconnected diagrams. It follows directlythat the connected contribution is also kernel-independent, because a conn µ = a hlbl µ − a disc µ inour SU(3) f case. One can further show that a conn µ is also itself kernel-independent in generalflavor cases. [1] G. W. Bennett et al. (Muon ( g − , 072003 (2006), arXiv:hep-ex/0602035.[2] F. Jegerlehner and A. Nyffeler, “The muon g − , 1 (2009), arXiv:0902.3360[hep-ph].[3] T. Blum, A. Denig, I. Logashenko, E. de Rafael, B. Lee Roberts, et al. , “The muon ( g − The Anomalous Magnetic Moment of the Muon (Springer Tracts Mod. Phys.,vol. 274, 2017).[5] T. Aoyama et al. , “The anomalous magnetic moment of the muon in the Standard Model,”(2020), arXiv:2006.04822 [hep-ph].[6] J. Grange et al. (E989), “Muon ( g −
2) technical design report,” (2015), arXiv:1501.06858[physics.ins-det].[7] T. Mibe (J-PARC g − g − Tau lepton physics. Proceedings, 11th International Workshop, TAU 2010,Manchester, UK, September 13-17, 2010 , Nucl. Phys. B (Proc. Suppl.) , 242 (2011). [8] H. B. Meyer and H. Wittig, “Lattice QCD and the anomalous magnetic moment of the muon,”Prog. Part. Nucl. Phys. , 46 (2019), arXiv:1807.09370 [hep-lat].[9] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, and C. Lehner, “Connectedand leading disconnected hadronic light-by-light contribution to the muon anomalous magneticmoment with a physical pion mass,” Phys. Rev. Lett. , 022005 (2017), arXiv:1610.04603[hep-lat].[10] N. Asmussen, A. G´erardin, A. Nyffeler, and H. B. Meyer, “Hadronic light-by-light scatter-ing in the anomalous magnetic moment of the muon,” SciPost Phys. Proc. , 031 (2019),arXiv:1811.08320 [hep-lat].[11] G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer, “Hadronic light-by-light contribu-tion to ( g − µ : a dispersive approach,” Proceedings, 35th International Symposium on LatticeField Theory (Lattice 2017): Granada, Spain, June 18-24, 2017 , EPJ Web Conf. , 01025(2018), arXiv:1711.00281 [hep-ph].[12] A. G´erardin, H. B. Meyer, and A. Nyffeler, “Lattice calculation of the pion transition formfactor with N f = 2 + 1 Wilson quarks,” Phys. Rev. D , 034520 (2019), arXiv:1903.09471[hep-lat].[13] J. Gronberg et al. (CLEO), “Measurements of the meson-photon transition form factors of lightpseudoscalar mesons at large momentum transfer,” Phys. Rev. D , 33 (1998), arXiv:hep-ex/9707031.[14] P. del Amo Sanchez et al. (BaBar), “Measurement of the γγ ∗ → η and γγ ∗ → η (cid:48) transitionform factors,” Phys. Rev. D , 052001 (2011), arXiv:1101.1142 [hep-ex].[15] J. P. Lees et al. (BaBar), “Measurement of the γ (cid:63) γ (cid:63) → η (cid:48) transition form factor,” Phys. Rev.D , 112002 (2018), arXiv:1808.08038 [hep-ex].[16] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, and C. Lehner, “Hadroniclight-by-light scattering contribution to the muon anomalous magnetic moment from latticeQCD,” Phys. Rev. Lett. , 132002 (2020), arXiv:1911.08123 [hep-lat].[17] N. Asmussen, A. G´erardin, J. R. Green, H. B. Meyer, and A. Nyffeler, “Hadronic light-by-light scattering contribution to the muon g − g − Proceedings, 37th International Symposium on Lattice Field Theory(Lattice 2019), Wuhan, China, 16–22 June 2019 , (2019), arXiv:1911.05573 [hep-lat].[19] A. G´erardin, T. Harris, and H. B. Meyer, “Nonperturbative renormalization and O ( a )-improvement of the nonsinglet vector current with N f = 2 + 1 Wilson fermions and tree-level Symanzik improved gauge action,” Phys. Rev. D , 014519 (2019), arXiv:1811.08209[hep-lat].[20] M. Bruno et al. , “Simulation of QCD with N f = 2 + 1 flavors of non-perturbatively improvedWilson fermions,” JHEP , 043 (2015), arXiv:1411.3982 [hep-lat].[21] M. Bruno, T. Korzec, and S. Schaefer, “Setting the scale for the CLS 2 + 1 flavor ensembles,”Phys. Rev. D , 074504 (2017), arXiv:1608.08900 [hep-lat].[22] M. Peardon, J. Bulava, J. Foley, C. Morningstar, J. Dudek, R. G. Edwards, B. Jo´o, H.-W. Lin, D. G. Richards, and K. J. Juge (Hadron Spectrum), “Novel quark-field creationoperator construction for hadronic physics in lattice QCD,” Phys. Rev. D , 054506 (2009),arXiv:0905.2160 [hep-lat]. [23] C. Morningstar, J. Bulava, J. Foley, K. J. Juge, D. Lenkner, M. Peardon, and C. H. Wong,“Improved stochastic estimation of quark propagation with Laplacian Heaviside smearing inlattice QCD,” Phys. Rev. D , 114505 (2011), arXiv:1104.3870 [hep-lat].[24] R. A. Briceno, J. J. Dudek, R. G. Edwards, and D. J. Wilson, “Isoscalar ππ scattering andthe σ meson resonance from QCD,” Phys. Rev. Lett. , 022002 (2017), arXiv:1607.05900[hep-ph].[25] A. Nyffeler, “Precision of a data-driven estimate of hadronic light-by-light scatteringin the muon g −
2: Pseudoscalar-pole contribution,” Phys. Rev.
D94 , 053006 (2016),arXiv:1602.03398 [hep-ph].[26] P. Masjuan and P. Sanchez-Puertas, “Pseudoscalar-pole contribution to the ( g µ − , 054026 (2017), arXiv:1701.05829 [hep-ph].[27] N. Asmussen, A. G´erardin, H. B. Meyer, and A. Nyffeler, “Exploratory studies for theposition-space approach to hadronic light-by-light scattering in the muon g − Proceed-ings, 35th International Symposium on Lattice Field Theory (Lattice 2017): Granada, Spain,June 18-24, 2017 , EPJ Web Conf. , 06023 (2018), arXiv:1711.02466 [hep-lat].[28] A. G´erardin, J. Green, O. Gryniuk, G. von Hippel, H. B. Meyer, V. Pascalutsa, and H. Wittig,“Hadronic light-by-light scattering amplitudes from lattice QCD versus dispersive sum rules,”Phys. Rev. D , 074501 (2018), arXiv:1712.00421 [hep-lat].[29] J. Bijnens, “Hadronic light-by-light contribution to a µ : extended Nambu-Jona-Lasinio, chiralquark models and chiral Lagrangians,” EPJ Web Conf. , 01002 (2016), arXiv:1510.05796[hep-ph].[30] J. Bijnens and J. Relefors, “Pion light-by-light contributions to the muon g − ,113 (2016), arXiv:1608.01454 [hep-ph].[31] J. Bijnens, E. Pallante, and J. Prades, “Analysis of the hadronic light-by-light contributionsto the muon g − , 379 (1996), arXiv:hep-ph/9511388.[32] M. Knecht, S. Narison, A. Rabemananjara, and D. Rabetiarivony, “Scalar meson con-tributions to a µ from hadronic light-by-light scattering,” Phys. Lett. B , 111 (2018),arXiv:1808.03848 [hep-ph].[33] J. Prades, E. de Rafael, and A. Vainshtein, “The hadronic light-by-light scattering contribu-tion to the muon and electron anomalous magnetic moments,” Adv. Ser. Direct. High EnergyPhys. , 303 (2009), arXiv:0901.0306 [hep-ph].[34] A. G´erardin, M. C`e, G. von Hippel, B. H¨orz, H. B. Meyer, D. Mohler, K. Ottnad, J. Wilhelm,and H. Wittig, “The leading hadronic contribution to ( g − µ from lattice QCD with N f = 2+1flavours of O( a ) improved Wilson quarks,” Phys. Rev. D , 014510 (2019), arXiv:1904.03120[hep-lat].[35] M. Hoferichter, B.-L. Hoid, B. Kubis, S. Leupold, and S. P. Schneider, “Pion-pole contributionto hadronic light-by-light scattering in the anomalous magnetic moment of the muon,” Phys.Rev. Lett. , 112002 (2018), arXiv:1805.01471 [hep-ph].[36] G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer, “Dispersion relation for hadroniclight-by-light scattering: two-pion contributions,” JHEP , 161 (2017), arXiv:1702.07347[hep-ph].[37] M. L¨uscher and S. Schaefer, “Lattice QCD with open boundary conditions and twisted-massreweighting,” Comput. Phys. Commun. , 519 (2013), arXiv:1206.2809 [hep-lat].[38] R. G. Edwards and B. Jo´o (SciDAC, LHPC, UKQCD), “The Chroma software system forlattice QCD,” Lattice field theory. Proceedings, 22nd International Symposium, Lattice 2004,Batavia, USA, June 21-26, 2004 , Nucl. Phys. B (Proc. Suppl.) , 832 (2005), arXiv:hep- lat/0409003.[39] A. Stathopoulos and J. R. McCombs, “PRIMME: PReconditioned Iterative MultiMethodEigensolver—methods and software description,” ACM Trans. Math. Softw. , 21:1 (2010).[40] J¨ulich Supercomputing Centre, “JUQUEEN: IBM Blue Gene/Q supercomputer system at theJ¨ulich Supercomputing Centre,” J. Large-Scale Res. Facil. , A1 (2015).[41] C. W. Bernard and M. F. Golterman, “Chiral perturbation theory for the quenched approxi-mation of QCD,” Phys. Rev. D , 853 (1992), arXiv:hep-lat/9204007.[42] S. R. Sharpe and N. Shoresh, “Physical results from unphysical simulations,” Phys. Rev. D , 094503 (2000), arXiv:hep-lat/0006017.[43] L. Giusti and M. L¨uscher, “Chiral symmetry breaking and the Banks-Casher relation in latticeQCD with Wilson quarks,” JHEP , 013 (2009), arXiv:0812.3638 [hep-lat].[44] M. Della Morte and A. J¨uttner, “Quark disconnected diagrams in chiral perturbation theory,”JHEP , 154 (2010), arXiv:1009.3783 [hep-lat].[45] J. Wess and B. Zumino, “Consequences of anomalous Ward identities,” Phys. Lett. B , 95(1971).[46] E. Witten, “Global aspects of current algebra,” Nucl. Phys. B , 422 (1983).[47] W. Detmold, B. Tiburzi, and A. Walker-Loud, “Electromagnetic and spin polarizabilities inlattice QCD,” Phys. Rev. D , 114505 (2006), arXiv:hep-lat/0603026.[48] S. Laporta and E. Remiddi, “The analytical value of the electron light-light graphs contribu-tion to the muon ( g −
2) in QED,” Phys. Lett. B301