Hadronic light-by-light contribution to the muon anomaly from lattice QCD with

infinite volume QED at physical pion mass

Thomas Blum,1, ∗ Norman Christ,2 Masashi Hayakawa,3, 4 Taku Izubuchi,5, 6
Luchang Jin,1, 6, † Chulwoo Jung,5 Christoph Lehner,7 and Cheng Tu1
(RBC and UKQCD Collaborations)
1
Physics Department, University of Connecticut, Storrs, CT 06269-3046, USA
2
Physics Department, Columbia University, New York, NY 10027, USA
3
Department of Physics, Nagoya University, Nagoya 464-8602, Japan
4
Nishina Center, RIKEN, Wako, Saitama 351-0198, Japan
5
Physics Department, Brookhaven National Laboratory, Upton, NY 11973, USA
6
RIKEN-BNL Research Center, Brookhaven National Laboratory, Upton, NY 11973, USA
7
Fakultät für Physik, Universität Regensburg, Universitätsstraße 31, 93040 Regensburg, Germany
arXiv:2304.04423v2 [hep-lat] 6 Dec 2024

(Dated: December 9, 2024)

The hadronic light-by-light scattering contribution to the muon anomalous magnetic moment,
(g − 2)/2, is computed in the infinite volume QED framework with lattice QCD. We report aHLbL µ =
12.47(1.15)(0.95) × 10−10 where the first error is statistical and the second systematic. The result is
mainly based on the 2+1 flavor Möbius domain wall fermion ensemble with inverse lattice spacing
a−1 = 1.73 GeV, lattice size L = 5.5 fm, and mπ = 139 MeV, generated by the RBC-UKQCD
collaborations. The leading systematic error of this result comes from the lattice discretization.
This result is consistent with previous determinations.

I. INTRODUCTION q = p′ − p, ν q = p′ − p, ν

Muons are spin-1/2 charged particles with non-zero

p′ p p′
magnetic moment: p

e
µ = −g S, (1) FIG. 1. The HVP (left) and HLbL (right) contribution to
2m muon g − 2. The shaded circles represent hadronic interac-
where S is the particle’s spin, e and m are the elec- tions.
tric charge and mass, respectively, and g is the Landé
g-factor. The Dirac equation predicts that g = 2, ex-
actly, so any difference from 2 must arise from interac- The muon anomalous magnetic moment is one of the
tions. The magnetic moment of a fermion can be defined most precisely determined quantities in particle physics.
in terms of its electromagnetic form factors in the zero Compared with the electron anomalous magnetic mo-
momentum transfer limit. Lorentz and gauge symme- ment, which is determined with higher accuracy, the
tries tightly constrain the form of the interactions. In muon is expected to be much more sensitive to new
Euclidean space-time: physics at very large energy scales due to its heavier
mass. Two experiments Fermilab (E989) [1] and J-PARC
hµ(p′ , s′ )|Jν (0)|µ(p, s)i (2) (E34) [2, 3] are aiming at even higher precision in mea-
 2
 suring the muon g − 2. The initial result released by
F2 (q )
= −eūs′ (p′ ) F1 (q 2 )γν + i [γν , γρ ]qρ us (p), Fermilab (E989) [1] confirmed the previously best result
4m
obtained by the BNL E821 experiment [4] and reduced
where Jν is the electromagnetic current, and F1 and F2 the experimental uncertainty from 0.54 ppm to 0.46 ppm.
are form factors, giving the charge and magnetic moment The final goal of the Fermilab experiment is to reduce the
at zero momentum transfer (q 2 = (p′ − p)2 = 0), or uncertainty further to approximately 0.14 ppm. The J-
static limit. The us (p) and ūs (p) are Dirac spinors. The PARC experiment adopts a very different measurement
anomalous part of the magnetic moment is given by F2 (0) strategy and will serve as an important cross-check. Its
alone, and is known as the anomaly, final accuracy goal is 0.45 ppm for statistical uncertainty
and 0.07 ppm for systematic uncertainty.
g−2 The Standard Model result provided by the Muon g −2
aµ ≡ F2 (q 2 = 0) = . (3)
2 Theory Initiative [5–25] currently has an uncertainty of
0.37 ppm and is in 4.2σ tension with the experimental
value. The two leading sources of uncertainty come from
the leading QCD contributions, hadronic vacuum polar-
∗ [email protected] ization (HVP) and hadronic light-by-light (HLbL) scat-
† [email protected] tering, both illustrated in Fig. 1.
 2

For a long time, the HLbL contribution has been esti- carefully calculated and were found to be consistent with
mated by hadronic models [26–29]. In Ref. [30, 31], the zero. The charm quark contribution from the connected
hadronic light-by-light scattering amplitudes estimated and disconnected diagrams were calculated as well. Fi-
by hadronic models are compared with lattice QCD cal- nally, a statistically more precise result was obtained:
culations. The result from hadronic models is difficult
to improve further. More recently, the HLbL contribu- aHLbL,Mainz2021
µ × 1010 = 10.96(1.59). (6)
tion has been calculated with the dispersion relation ap-
In this work, we present our latest lattice calculation of
proach [17–25, 32–37]. These results are summarized in
the HLbL contribution to muon g−2 with subtracted infi-
the 2020 white paper [5],
nite volume QED weighting function as developed in our
aHLbL,WP2020,pheno × 1010 = 9.2(1.9) (4) previous work [46]. The calculation is directly performed
µ
at physical pion mass and therefore eliminates the sys-
Within the dispersion relation framework, the pseudo- tematic uncertainty from chiral extrapolations. We find
scalar pole is the leading contribution to HLbL. This that calculating the HLbL contribution directly at the
contribution is defined in terms of the pseudo-scalar tran- physical pion mass is considerably more difficult than a
sition form factors to two (possibly off-shell) photons, calculation at a heavier pion mass due the larger contri-
π 0 , η, η ′ → γ ∗ γ ∗ . These form factors can also be calcu- bution and statistical fluctuations from the long distance
lated with lattice QCD. [21, 38–40] region. Following the suggestion of Ref. [48], we focus
The first direct lattice QCD calculation of the HLbL on the calculation of the connected and leading discon-
contribution was performed by the RBC-UKQCD col- nected diagrams. The main calculation is performed with
laborations. [41] Then, we made many improvements to a single lattice spacing where a−1 = 1.73 GeV. While we
the calculation method and also included the contribu- use Möbius domain wall fermions (MDWF) to eliminate
tion from the leading disconnected diagrams. [42, 43] In O(a) discretization errors, the remaining O(a2 ) effects
our previous work [24], we applied the method developed are estimated to be the largest source of systematic un-
to lattice gauge ensembles with different lattice spacings certainty of this calculation.
and spatial volumes, all at physical pion mass. The cal- The paper is organized as follows: in section II and
culations incorporate both QCD and QED in the same III, we review the calculation method developed in our
finite volume lattice using the QEDL scheme [44]. We previous works and describe new techniques used in this
extrapolated the results to infinite volume and the con- calculation. We then show our numerical results in sec-
tinuum limit. The final result was tion IV and conclude in section V.

aHLbL,RBC-UKQCD2019
µ × 1010 = 7.87(3.06)stat(1.77)syst .
(5) II. THEORETICAL FRAMEWORK AND
CALCULATION METHOD
This was the first lattice QCD result for HLbL with all
systematic effects controlled. The statistical error was We base our current calculation on the framework al-
still larger than the phenomenological results obtained ready set up in our previous works [42, 43, 46]. In this
with dispersion relations or the previous hadronic model work, we combine the method used in Ref. [42, 43] to
estimation. calculate the hadronic four point function with lattice
The Mainz group pioneered the QED∞ approach and QCD with the subtracted infinite volume QED weight-
pre-calculated the QED kernel semi-analytically in in- ing function developed in Ref. [46]. We have tested the
finite volume [30, 45]. We also explored this strategy QED weighting function by calculating the leptonic light-
and independently calculated the infinite volume QED by-light contribution to muon g − 2 and have reproduced
weighting function [46]. In that work, we found that the results obtained by analytical calculation [50, 51]. We
while the QED∞ approach ensures an exponentially sup- outline the method below and focus on the improvements
pressed finite volume error, the size of the finite volume we made in this work.
error can still be significant. However, making use of the All of the diagrams contributing at O(α3 ) to HLbL
current conservation property of the hadronic four point scattering are shown in Fig. 2. We use the term (quark-)
function, we designed subtractions to the infinite volume “connected diagram” to refer to the one on the left on
QED weighting functions to reduce the finite volume er- the top row and “disconnected diagram” to refer to the
rors and also the discretization errors. The subtracted one on the right. The remaining diagrams, which are all
infinite volume QED weighting function is used in the suppressed in the flavor SU(3) limit, are referred to as the
current work. “sub-leading disconnected diagrams”. We only explicitly
The effectiveness of the subtraction is confirmed by draw quark loops that are connected to photons. Gluons
the Mainz group and a different subtraction scheme was and sea quark loops that are not connected to photons are
adopted in their calculation. [45, 47–49]. Their calcu- not shown in the figure but are included automatically in
lation was performed with pion masses ranging from dynamical lattice QCD calculations. To make non-zero
200 MeV to 420 MeV and extrapolated to the physi- contributions to the muon g − 2, the quark loops in the
cal point. All sub-leading disconnected diagrams were disconnected diagrams must be connected by gluons.
 3

x y

tsrc α, ρ η, κ β, σ tsnk

FIG. 3. Diagramatic representation of the QED weighting

function defined in Eq. 9, following Ref. [46].

As is well known, the above expression contains an in-

FIG. 2. Diagrams contributing to the muon anomaly. The di- frared divergence that vanishes after projection to its
agrams in the top row are the leading ones, and do not vanish magnetic part. We can also remove this infrared diver-
in the SU (3) flavor limit. Strong interactions to all orders, gent piece by the following procedure:
including gluons connecting the quark loops and sea quark
loops which are not connected by photons, are not shown. 1
G(1)
σ,κ,ρ (y, z, x) = Gσ,κ,ρ (y, z, x) (11)
2
1
The contribution to the muon g − 2 can be calculated + [Gρ,κ,σ (x, z, y)]† .
2
with the combination of the hadronic four-point function
H and the QED weighting function G [46]: In addition we can perform somewhat arbitrary subtrac-
e Σi tions to this infinite volume QED weighting function
aHLbL
µ ūs′ (~0) us (~0) (7) without changing the final result due to vector current
m 2
conservation satisfied by the hadronic four point func-
1 XX 1

= ǫi,j,k xop − xref (x, y, z) j tion,
V T x x,y,z 2
op

× i e Hk,ρ,σ,λ (xop , x, y, z)ūs′ (~0)Gρ,σ,λ (x, y, z)us (~0),

3 6 G(2) (1)
σ,κ,ρ (y, z, x) = Gσ,κ,ρ (y, z, x) (12)
− G(1)
σ,κ,ρ (z, z, x) − G(1)
σ,κ,ρ (y, z, z).
where ūs′ (~0), us (~0) are Dirac spinors for the outgoing and
incoming muon in the diagram, respectively. V stand for (1)
the spatial volume of the lattice and T stand for the size Note that Gσ,κ,ρ (z, z, z) = 0, so this subtraction signif-
of the temporal extent of the lattice. Σk = ǫi,j,k γi γj /(2i) icantly reduces the size of the QED weighting function
is the 4 × 4 version of the Pauli matrix, σk . From the when |x − z| or |y − z| is small. This is the region where
spin structure of the muon particle, we can obtain the the hadronic function from the lattice calculation has the
expression for aHLbL : largest discretization error. It turns out that this subtrac-
µ
tion greatly reduces the discretization error. This is the
2me2 1 X X 1
major finding of Ref. [46]. We should also note that the
aHLbL
µ = ǫi,j,k xop − xref (x, y, z) j subtraction does impact the integrand and partial sum.
3 V T x x,y,z 2
op In Ref. [48], a modified subtraction scheme is used, so
× (6e4 )Hk,ρ,σ,λ (xop , x, y, z)Mi,ρ,σ,λ (x, y, z) (8) their integrand and partial sum cannot be directly com-
pared with ours. Finally, we include all possible permu-
where tations of the subtracted QED weighting function which
1 h1 3 i are required for the total contribution to the muon g − 2:
Mi,ρ,σ,λ (x, y, z) = Tr i Gρ,σ,λ (x, y, z)Σi . (9)
2 6
i3 Gρ,σ,κ (x, y, z) = G(2) (2)
ρ,σ,κ (x, y, z) + Gσ,κ,ρ (y, z, x) (13)
The QED weighting function G is shown diagramatically
in Fig. 3 and is expressed in terms of the free muon and + G(2) (2)
κ,ρ,σ (z, x, y) + Gκ,σ,ρ (z, y, x)
Feynman gauge photon propagators, Sµ and G: + G(2) (2)
ρ,κ,σ (x, z, y) + Gσ,ρ,κ (y, x, z).

Gσ,κ,ρ (y, z, x) = lim emµ (tsnk −tsrc ) (10)

tsrc →−∞,tsnk →∞ Another component of the master formula Eq. 8 is xref ,
Z
the reference position for the moment method to calcu-
× G(x, α)G(y, β)G(z, η) late the magnetic moment. Again, due to current con-
α,β,η,~
xsnk ,~
xsrc
servation, there are many possible choices for xref . In
× Sµ (xsnk , β) iγσ Sµ (β, η)iγκ Sµ (η, α)iγρ Sµ (α, xsrc ) . this work, we use the following choice for the connected
 4

diagram: 6e4 Hν,ρ,σ,κ

discon
(xop , x, y, z)
⇒ 6e4 Hν,ρ,σ,κ
discon-no-perm
(xop , x, y, z) (19)
xref (x, y, z) = xref-far (x, y, z) (14) D X  
 =3 e2q′ Tr γν Sq′ (xop , x)γρ Sq′ (x, xop )
 x if |y − z| < min(|x − y|, |x − z|) q′ =u,d,s
 
y if |x − z| < min(|x − y|, |y − z|) X
= × e2q Tr γκ Sq (z, y)γσ Sq (y, z)
 z1
 if |x − y| < min(|x − z|, |y − z|)
(x + y + z) otherwise q=u,d,s
3 E
− γκ Sq (z, y)γσ Sq (y, z) QCD
,
QCD
We make a slightly different choice of xref , Eq. (21), for
the disconnected diagrams. The rationale will be de- where Sq denotes the quark propagator. As shown in
scribed in the later part of this section. We use H to Fig. 4, we calculate the hadronic four-point function with
denote the hadronic four point function: two point-source propagators. The above trick can make
the evaluation more efficient [42, 47]. We note that
(6e4 )Hk,ρ,σ,λ (xop , x, y, z)
= hT Jk (xop )Jρ (x)Jσ (y)Jλ (z)iQCD (15) xop , ν

X z, κ
Jν (x) = eq ZV ψ̄q (x)γν ψq (x) (16) y, σ x, ρ
q=u,d,s,c

where ZV is the lattice local vector current renormal-

ization constant. After Wick contraction, H can be ex- xsrc xsnk
pressed as the sum of different types of diagrams as il-
lustrated in Fig. 2. For the disconnected and sub-leading xop , ν
disconnected diagrams, we require the quark loops to be
connected by gluons.
x, ρ y, σ z, κ
con discon sub-leading-discon charm
H=H +H +H +H (17)

where Hcon includes the light and strange quark con-

nected diagrams, Hdiscon the light and strange quark xsrc xsnk
disconnected diagrams, Hsub-leading-discon the light and
strange quark sub-leading disconnected diagrams (van- FIG. 4. Two point-source propagators are used to calculate
ish in the flavor SU(3) limit), and Hcharm all diagrams the hadronic four-point function. The locations of the point-
involving charm quark loops. sources are indicated in the diagrams as small circles.
Naturally, after the Wick contraction, H defined in
Eq. (15) can be expressed in terms of quark propagators Hdiscon-no-perm, the disconnected diagram of the hadronic
and includes all permutations of x, y, and z. However, four-point function without permutation, still satisfies
note that all other terms in the master formula Eq. (8) the current conservation condition in the continuum limit
are symmetric with respect to permutations of x, y, and (this is not true for Hcon-no-perm ) which leads to the fol-
z (along with its Lorentz indices). Therefore, we are lowing relation in the infinite volume limit:
allowed to calculate only a subset of the diagrams gener- X discon-no-perm
ated by the Wick contractions in Eq. (15) and multiply Hk,ρ,σ,λ (xop , x, y, z) = 0. (20)
them with appropriate factors. xop

We are therefore allowed to alter the choice of xref for the

6e4 Hν,ρ,σ,κ
con
(xop , x, y, z) disconnected diagram:
⇒ 6e4 Hν,ρ,σ,κ
con-no-perm
(xop , x, y, z) (18) xref-discon = x. (21)
D X 
4
= −6 Re eq Tr γν Sq (xop , x)γρ Sq (x, z) This choice allows the summation over xop to be per-
q=u,d,s formed independently of coordinates y and z. This
E
× γκ Sq (z, y)γσ Sq (y, xop ) , choice also has the benefit of suppressing the contribution
QCD from the region where |x − xop | is small and the quark
 5

loop is large before subtracting the vacuum expectation speed up the sloppy propagator calculation, we use the
value. Also, we note that the new choice xref-discon is zMöbius domain wall fermion formulation [54] with re-
the same as the initial choice in the long distance region duced fifth dimension size Ls = 14 to approximate the
where the contribution is dominated by π 0 -exchange and unitary MDWF action used in the gauge ensemble gen-
|y − z| ≪ min(|x − y|, |x − z|). This property guaran- eration (Ls = 24). For deflation of the Dirac operator
tees that the connected and disconnected diagrams have we reuse the eigenvectors generated for the calculation
exactly the same QED weighting function in the long of the hadronic vacuum polarization contribution to the
distance region, which leads to Eq. 28 before taking the muon g −2 [55, 56] using the locally-coherent Lanczos ap-
continuum limit. This is the main reason for the choice proach [57]. Then, we use a two-level all-mode-averaging
of xref for the connected diagrams. (AMA) method [58] to correct the bias caused by the
inexact propagators. The AMA method requires a por-
tion of the sloppy propagators be computed again with
III. LATTICE DETAILS higher precision. In the first step a more precise version
of the light and strange quark propagators are calculated
The calculation presented here is performed on ensem- on a randomly selected set of 64 points among the 1024
bles of gauge fields generated by the RBC and UKQCD points. In the second step among the 64 points, we ran-
collaborations [52]. The main calculation is carried out domly select 16 points and calculate these propagators
using the 48I ensemble, a physical-mass ensemble gen- to full precision. These propagators are used to calculate
erated with 2+1 flavors of Möbius domain wall fermions both the connected and the disconnected diagrams in the
(MDWF). We also use a few other ensembles to calculate bias correction. By using two steps, we are able to reduce
various corrections and estimate systematic effects. The the statistical error coming from the bias correction to a
relevant information about the 48I ensemble and other level well below the statistical error for the sloppy part.
ensembles is listed in Table I. We always use the MDWF For the connected diagrams, we sample two-point pairs
action. For most of the ensembles, the quarks have their among all the possible 2048 × (2048 − 1)/2 + 2048 =
physical masses. 2098176 combinations (including the case where x = y).
For each of the two-point pairs, we can perform the con-
48I 64I 24D 32D 24DH traction as described in Eq. (18) in the previous section.
mπ (MeV) 139 135 142 142 341 As the number of all possible pairs is enormous, the cost
a−1 (GeV) 1.730 2.359 1.015 1.015 1.015 to perform all the contractions is not affordable. There-
a (fm) 0.114 0.084 0.194 0.194 0.194 fore, we only calculate the contraction for a subset of
L (fm) 5.47 5.38 4.67 6.22 4.67 the available two-point pairs. We sample the subset with
Ls 24 12 24 24 24 the empirical probability p(r), which is a function of the
# meas 113 64 31 70 37 distance between the two points, r = |x − y|,
TABLE I. 2+1 flavors of MDWF gauge field ensembles gen-  Npsrc −1

 3 if r = 0
erated by the RBC-UKQCD collaborations [52]. The labels  2(L T −1)
indicate the lattice size in lattice units and the QCD gauge 1 if 8 ≥ r > 0
p(r) = 1 , (22)
action, where “I” stands for Iwasaki, and “D” stands for 
 if L ≥ r > 8
 (r/8)3
Iwasaki+DSDR. DSDR stands for dislocation-suppressing- 0 if r > L
determinant-ratio, which is an additional term in the action
to soften explicit chiral symmetry breaking effects, needed in where L = 48, T = 96, and the numbers are in lattice
particular for very coarse lattices [53]. The pion mass mπ , units. On average, about 57, 000 light quark point-source
lattice spacing aa , spatial extent L, extra fifth dimension size propagator pairs per configuration are sampled to calcu-
Ls , and the number of QCD configurations used are listed. late the connected diagrams. The reason we sample long-
distance point-pairs less often is because the connected
a After the analysis for this work was completed, we updated the hadronic four point function and its statistical fluctua-
lattice spacing for the 24D and 32D ensembles. The new lattice tion decrease with distance for long distances. For each
spacing is a−1 = 1.023 GeV and makes no material difference to
pair, due to the sampling procedure, the probability of
the results presented here.
that pair taking a particular relative coordinate is a func-
tion of r. The total connected diagram contribution is
For our main calculation on the 48I ensemble, we use equal to the expectation value of the contraction for the
113 configurations. On each configuration, we randomly point pairs multiplied by the inverse probabilities. We
sample 2048 uniformly distributed distinct points and refer to the inverse probability as the weight w(r):
calculate light quark point-source propagators for each

of the selected points. w0 if r = 0
Among the above 2048 points we randomly select w(r) = , (23)
w0 /p(r) if r > 0
1024 points and also calculate strange quark propaga-
tors. These 2048 + 1024 point-source propagators are X
computed using the (deflated-) conjugate gradient al- w0 = p(|x|). (24)
gorithm with a sloppy stopping condition. To further x
 6

To further reduce the contraction cost and the storage For the light quark loop, there are 9, 156, 431 points
cost of saving these propagators, we use field sparsening within the |z − y| ≤ L range. Among these, we sam-
techniques [59, 60]. We randomly sample 1/16 points ple about 17, 000 points for z on average. For each of the
among all the points of the lattice and only perform con- sampled z, we then loop over all possible x locations with
tractions on these randomly selected points. Since we |x − y| ≤ L and perform the full contraction. Due to the
only need propagator values on these points, we only save sampling procedure, we multiply the result of the con-
these values to disk. To accommodate the sparsening, we traction by the inverse of the sample probability 1/py (z)
multiply a factor of 162 for the summation over z and xop to obtain an unbiased final result. While the sampled z
in Eq. (8), except when z = xop , where we multiply by 16. points amount to less than 0.2% of the possible z points
In this way the sampling procedure does not introduce in the allowed range, we expect the final statistical pre-
any systematic effects on our central value. cision to be almost unaffected by the adaptive sampling
For the disconnected diagrams, the two point-source procedure. That is, we expect the final statistical error
locations are also indicated in Fig. 4. Different from would be almost the same even if we had calculated the
the connected diagrams, the statistical fluctuations of the contraction using all the z locations. This expectation
disconnected hadronic functions are not suppressed when is based on a quick test run with much lower threshold
r = |x − y| increases. Therefore, we use all possible com- t0 = 10−3 , which corresponds to about a factor of 20
binations of point pairs of x and y as long as |x−y| ≤ L to smaller number of sampled points; we find the final sta-
estimate the summation over x and y in Eq. (8). Thanks tistical error is roughly the same.
to the choice of xref = xref-discon = x in Eq. (21), the sum- In contrast to the connected diagram, we do not
mation over xop can be performed for each point-source sparsen the propagators. To save disk space, we tem-
propagator with source location at x, the result can be porarily save the following intermediate contraction for
used for all possible values of y locations. However, the each point-source propagator:
summation of z depends on the location of x and y due  
to the QED weighting function. Therefore, we need to Tr γκ Sq (z, y)γσ Sq (y, z) . (26)
perform the summation over z for each pair of points.
To make the contractions of all point pairs affordable, For each propagator on each lattice site, the above con-
we aggressively sparsen when summing over z. Fortu- traction corresponds to 4 × 4 complex numbers while the
nately, the vertex z and vertex y are connected by two full propagator takes 12 × 12 complex numbers.
quark propagators, and the four point function is expo-
nentially suppressed when the distance between z and y
increases. IV. RESULTS
Naturally, we can sample the z locations based on the
distance between z and y, similar to the sampling of x In this section we display several figures of the sum-
and y combinations in the calculation of the connected di- mand (integrand) or its partial-sum as a function of
agrams. However, for this disconnected diagram, we dis- Rmax :
covered a much more efficient adaptive sampling scheme
as follows. For each point-source location y, we calculate Rmax = max(|x − y|, |y − z|, |x − z|), (27)
the following square norm n(z, y) of the quark loop for
all z: where x, y, and z represent the positions of the three in-
X  ternal quark-photon vertices. Partial-sums are obtained
n(z, y) = Tr γκ Sq (z, y)γσ Sq (y, z) by summing the corresponding summand from 0 up to
κ,σ the specified value of Rmax (inclusive). The rightmost
 2 value in a partial-sum plot corresponds to the total con-
− γκ Sq (z, y)γσ Sq (y, z) QCD tribution.
Since the summand decreases exponentially at large
Unlike the full contraction, the above norm is indepen- Rmax , we expect the partial-sum to approach a plateau
dent of the location of x, and the location of z is sampled for large enough Rmax . We usually use the result of the
with the following (empirical) probability partial-sum at 4 fm, which includes all contributions from

 1p if n(z, y) ≥ t20 and |z − y| ≤ L the region Rmax < 4 fm. In the lattice calculation, we
py (z) = n(z, y)/t0 if n(z, y) < t20 and |z − y| ≤ L , calculate the partial sum with respect to Rmax and save
 the partial sum for all integer values (in lattice units)
0 if |z − y| > L
of the upper limit of Rmax . To compare with results
(25)
calculated using ensembles with different lattice spacings,
where t0 is the threshold parameter to control the overall we can linearly interpolate the partial sum for any upper
sample frequency. With this probability distribution, we limit of Rmax . In this work, we choose to interpolate at
sample more often points where the (norm of) the quark integer multiples × 0.1 fm. The summand is obtained at
loop is large which is much more efficient than basing half-integer multiples × 0.1 fm by taking the difference of
the probability on the distance |z − y|. We choose t0 = the linearly interpolated partial sum. Therefore, the data
5 × 10−5 in this work for the 48I ensemble. points of the plots shown in this work always have spacing
 7

0.1 fm, and the data points of the integrand plots always 48I light no-pion summand
fit1
have its x-axis value equal to half-integer times 0.1 fm. 15
In this work, we present our results for the anomalous
magnetic moment in the unit of 10−10 .
10

A. Light quark contribution

aµ × 1010 summand
5

48I light con summand

48I light discon summand 0
48I light summand
15

−5
10
aµ × 1010 summand

5 −10
0 1 2 3 4 5
Rmax (fm)
0
48I light no-pion
fit1
12
−5

10
−10

8
−15 aµ × 1010 partial-sum
0 1 2 3 4 5
6
Rmax (fm)
48I light con
48I light discon 4
48I light
25
2

20
0

15
aµ × 1010 partial-sum

−2
0 1 2 3 4 5
10 Rmax (fm)

5
FIG. 6. Light quark contributions from a special combina-
tion of the connected and disconnected diagrams to cancel
0 the long-distance π 0 exchange contribution (see Eq. 29). The
upper plot shows the summand and the lower plot shows the
−5 partial sum. The combined summand vanishes much faster
than the light quark diagrams by themselves. Solid curves
−10 represent a fit of the data to Eq. (30). The fit starts at 0.5 fm.
0 1 2 3 4 5 The dashed lines denote the statistical uncertainty of the fit.
Rmax (fm)

FIG. 5. Light quark contributions computed on the 48I en-

semble from the connected diagrams, the disconnected dia- gions where Rmax is as large as 3 or 4 fm. The statistical
grams, and the total. The upper plot shows the correspond- error is quite significant. Making the situation worse,
ing summands and the lower plot shows the partial sum. the contributions from the long-distance region of the
connected and the disconnected diagrams are opposite
The contributions from the connected diagrams, the in sign, while the statistical error of the connected and
disconnected diagrams, and the sum of the two contribu- the disconnected diagrams are largely independent. As a
tions are shown in Fig. 5. As can be seen from the figure, result, the relative uncertainty on the sum is much larger
the statistical error grows as Rmax increases. However, than that for the individual contributions. Results are
for both the connected and the disconnected diagrams, given in Tab. II as “48I light con Rmax < 4fm”, “48I
there are non-negligible contributions that come from re- light discon Rmax < 4fm”, and “48I light Rmax < 4fm”.
 8

To improve the situation, note that the large contri- The advantage comes in using the early plateau value
bution and the cancellation between the connected and of ano-pion
µ and combining it with acon µ , which plateaus
disconnected diagrams at long distance are due to π 0 ex- at much larger Rmax . This is a “hybrid” approach to
change and are well-understood theoretically [48, 61, 62]. calculate adiscon
µ and atotal
µ , which has a much smaller
It has been shown that, at large Rmax , the ratio of the statistical error than the direct combination. For acon µ ,
disconnected and the connected hadronic four-point func- we will use 4 fm as the upper limit of Rmax . The results
tion is −25/34. In our present computational setup, we are shown in Tab. II as entries that start with “48I light
use the same infinite volume QCD weighting function for discon Rmax < 4fm hybrid-” and “48I light Rmax < 4fm
both the connected and disconnected diagrams, and use hybrid-”.
the same variable Rmax to study the partial-sum of the The fitting function form in Eq. (30) is inspired by the
connected and disconnected diagrams.[63] Therefore the function in Eq. (20) of Ref. [48]. The functional form
same ratio applies to the contribution to aµ . Formally, used here differs due to the meanings of Rmax and the
we have: variable “|y|” in Ref. [48]. Also, note that the subtrac-
adiscon (Rmax > R) 1 (e2 + e2d )2 25 tion scheme of the QED weighting function is also differ-
µ
lim = − · u4 =− . ent [46–48, 65]. We tried using this fit function to fit other
con
R→∞ aµ (Rmax > R) 2 eu + e4d 34 contributions. The results are shown in Appendix C.
(28) The statistical errors are estimated by assuming the
Here, Rmax is defined in Eq. (27), being a function of results from each 48I configuration that we performed
the three vertex locations x, y, z as shown in Fig. 4. the measurements are statistically independent. We also
The long-distance contributions for the connected and tried to bin the data with different lengths in MD time
disconnected diagrams are denoted as acon units. The results are shown in Table. III. We can see
µ (Rmax > R)
and aµdiscon
(Rmax > R). They can be calculated with the that the statistical error does not significantly depend on
same summation as in Eq. (8) but with the additional the binning sizes. More interestingly, the statistical error
constraint that Rmax > R, where R is the long-distance calculated by assuming all the individual samples in each
configuration are independent is only a little bit smaller
cutoff. Note we take the R → ∞ limit in the above equa-
tion. This ratio is exact in this limit and is not affected than the statistical error calculated from the fluctuation
of different configurations. This implies that the current
by lattice artifacts or finite volume effects. Therefore, we
can construct the following combination: statistical errors still mainly come from the sampling of
the points from each configuration, instead of from the
25 con gauge field fluctuations sampled by different gauge con-
ano-pion
µ = adiscon
µ + a , (29)
34 µ figurations. Also, this shows that the strategy using dif-
where the π 0 exchange contribution to ano-pion vanishes ferent combinations of the point source propagators as
µ
in the long distance. We plot the summand and partial- point pairs to calculate the connected diagrams is very
sum of ano-pion in Fig. 6. Indeed, the partial sum of effective.
µ
ano-pion reaches the plateau much earlier than the con- The remaining finite volume, pion mass, and non-
µ
nected or disconnected diagrams. This trick of combining zero lattice spacing corrections will be studied in Sec-
the connected and disconnected with appropriate factors tions IV D, IV E, and IV F.
to obtain a faster plateau was employed in Ref. [64]. For
ano-pion
µ , we will use 2.0 fm and 2.5 fm as the upper limit
of Rmax . The results are shown in Tab. II as “48I light no- B. Strange quark contribution
pion Rmax < 2.5fm” or “48I light no-pion Rmax < 2fm”.
In the upper panel of Fig. 6, we fit the summand to The contributions from the strange quark connected
the following empirical form: diagrams, the disconnected diagrams, and the sum of the
6 two contributions are plotted in Fig. 7. Note the strange
Rmax
f (Rmax ) = A/fm4 3
e−BRmax /(fm·GeV) quark disconnected diagrams include diagrams where one
Rmax + (C fm)3 or both loops are strange quark loops, while the light
(30) quark disconnected diagrams discussed in the previous
The fit range is from 0.5 fm to 4 fm. We use the result section contain only light quarks. As seen in the fig-
of the fit to estimate the long-distance contribution to ure, the contribution from the strange quark-connected
ano-pion . Since this is a completely empirical fit, we will diagrams is very precise and very small. It can also be
µ
assign a 100% systematic uncertainty to it. The “no- clearly seen that the strange quark contribution vanishes
pion" results are also collected in Tab. II. much faster at long-distance compared to the light quark
With the definition of ano-pion , we obtain: contribution.
µ
The disconnected diagrams, on the other hand, are
25 con much noisier. However, we still expect the signal to van-
adiscon
µ = ano-pion
µ − a , (31)
34 µ ish faster than the light quark contribution due to the ab-
9 sence of the π 0 exchange contribution. Therefore, we can
atotal
µ = ano-pion
µ + acon . (32)
34 µ treat the strange quark disconnected diagrams as similar
 9

Contribution name aµ ×1010

48I light con Rmax < 2fm 7.24 (0.15)stat
48I light con Rmax < 2.5fm 10.64 (0.31)stat
48I light con Rmax < 4fm 18.61 (1.22)stat
48I light con Rmax > 4fm 7.56 (0.42)stat (1.06)syst [1.14]
48I light con FV-corr −1.79 (0.42)syst
48I light con mπ -corr 1.32 (0.27)stat (0.66)syst [0.72]
48I light con a2 -corr 0.00 (1.49)syst
light con 25.70 (1.33)stat (1.99)syst [2.39]
48I light no-pion Rmax < 2.5fm 5.09 (0.76)stat
48I light no-pion Rmax < 2fm 4.65 (0.42)stat
48I light no-pion Rmax > 2.5fm 0.31 (0.22)stat (0.31)syst [0.38]
48I light no-pion Rmax > 2fm 0.88 (0.46)stat (0.88)syst [0.99]
48I light discon Rmax < 2fm −0.67 (0.41)stat
48I light discon Rmax < 2.5fm −2.73 (0.74)stat
48I light discon Rmax < 4fm −7.49 (1.82)stat
48I light discon Rmax < 4fm hybrid-2.5fm −8.28 (1.31)stat (0.31)syst [1.35]
48I light discon Rmax < 4fm hybrid-2fm −8.15 (1.24)stat (0.88)syst [1.51]
48I light discon Rmax > 4fm −5.56 (0.31)stat (0.78)syst [0.84]
48I light discon FV-corr 1.31 (0.31)syst
48I light discon mπ -corr −0.98 (0.20)stat (0.49)syst [0.53]
48I light discon a2 -corr 0.00 (0.66)syst
light discon −12.71 (1.87)stat (1.17)syst [2.20]
light discon hybrid-2.5fm −13.50 (1.36)stat (1.21)syst [1.82]
light discon hybrid-2fm −13.37 (1.29)stat (1.46)syst [1.95]
48I light Rmax < 2fm 6.57 (0.43)stat
48I light Rmax < 2.5fm 7.90 (0.78)stat
48I light Rmax < 4fm 11.11 (2.11)stat
48I light Rmax < 4fm hybrid-2.5fm 10.32 (0.99)stat (0.31)syst [1.04]
48I light Rmax < 4fm hybrid-2fm 10.46 (0.89)stat (0.88)syst [1.25]
48I light Rmax > 4fm 2.00 (0.11)stat (0.28)syst [0.30]
48I light FV-corr −0.47 (0.11)syst
48I light mπ -corr 0.35 (0.07)stat (0.17)syst [0.19]
48I light a2 -corr 0.00 (0.83)syst
light total 12.99 (2.11)stat (0.90)syst [2.29]
light total hybrid-2.5fm 12.20 (1.01)stat (0.95)syst [1.38]
light total hybrid-2fm 12.33 (0.90)stat (1.25)syst [1.55]

TABLE II. Light quark contributions computed on the 48I ensemble. Values are different contributions to aµ × 1010 where
aµ = (gµ − 2)/2. The numbers in the square bracket are the statistical and systematic uncertainty combined in quadrature.
The tag “hybrid-” indicates the same quantity obtained using the hybrid approach.

Contribution name aµ ×1010 Bin(20) Bin(30) Bin(40) Bin(sample)

48I light con Rmax < 4fm 18.61 (1.22) 1.30 1.15 1.24 1.15
48I light discon Rmax < 2.5fm -2.73 (0.74) 0.79 0.84 0.75 0.65

TABLE III. Two light quark contributions computed on the 48I ensemble. The table lists the statistical errors calculated with
different binning sizes in MD time units in the HMC evolutions. The separation between the configurations we performed the
measurements is 10 MD time units. The last column named “Bin(sample)” list the statistical error calculated by assuming
all the individual samples in each configuration are independent. The samples for the connected diagrams are the point pairs,
where we have about 57,000 pairs for each configuration. The samples for the disconnected diagrams are the points for y in
Fig. 4 (summation over x is already performed), where we have about 2048 points for each configuration.

to ano-pion
µ . We listed the results with 4.0 fm, 2.5 fm, and con Rmax >”. The size of the long-distance connected
2 fm as the upper limit of Rmax . These results are shown contribution is used as the estimation for the systematic
in Tab. IV as entries that start with “48I strange discon uncertainty for the corresponding disconnected contribu-
Rmax <”. We estimate the systematic error caused by tions.
truncating the strange quark disconnected diagram in-
tegration using the size of the long-distance contribution
from the strange quark connected diagrams, which is also The remaining finite volume, pion mass, and non-
shown in Tab. IV as entries that start with “48I strange zero lattice spacing corrections will be studied in Sec-
tion. IV D, IV E, and IV F.
 10

Contribution name aµ ×1010

48I strange con Rmax < 4fm 0.327 (0.002)stat
48I strange con a2 -corr 0.026 (0.008)stat
strange con 0.353 (0.007)stat
48I strange con Rmax > 2.5fm 0.006 (0.000)stat
48I strange con Rmax > 2fm 0.024 (0.000)stat
48I strange discon Rmax < 4fm −0.380 (0.607)stat
48I strange discon Rmax < 2.5fm −0.357 (0.223)stat
48I strange discon Rmax < 2fm −0.280 (0.121)stat
48I strange discon a2 -corr 0.000 (0.029)syst
strange discon −0.380 (0.607)stat (0.029)syst [0.608]
strange discon hybrid-2.5fm −0.357 (0.223)stat (0.029)syst [0.225]
strange discon hybrid-2fm −0.280 (0.121)stat (0.037)syst [0.126]
48I strange Rmax < 4fm −0.053 (0.607)stat
48I strange Rmax < 4fm hybrid-2.5fm −0.030 (0.222)stat (0.006)syst [0.223]
48I strange Rmax < 4fm hybrid-2fm 0.047 (0.121)stat (0.024)syst [0.123]
48I strange a2 -corr 0.026 (0.008)stat (0.029)syst [0.030]
strange total −0.027 (0.607)stat (0.029)syst [0.608]
strange total hybrid-2.5fm −0.004 (0.223)stat (0.029)syst [0.225]
strange total hybrid-2fm 0.073 (0.121)stat (0.037)syst [0.127]

TABLE IV. Strange quark contributions computed on the 48I ensemble. Values are different contributions to aµ × 1010 where
aµ = (gµ − 2)/2. The numbers in the square bracket are the statistical and systematic uncertainty combined in quadrature.
The “strange discon” contribution includes disconnected diagrams where one or two loops are strange quark loops. The tag
“hybrid-” indicates the same quantity obtained using the hybrid approach.

C. Long-distance π 0 -exchange contribution The normalization constant Zπ is determined by the fol-

lowing requirement:
The dominant source of the HLbL contribution is ex-
pected to be the so-called “π 0 -pole” contribution [5]. Due h0|π 0 (x)|π 0 (~
p)i = eip·x . (37)
to the small mass of the π 0 , this contribution can be non-
negligible for large Rmax in a Euclidean space-time lattice We can then rewrite Eq. (33) in a slightly different form:
calculation, even where long-distance contributions and
finite volume effects are exponentially suppressed since it h0|T Jµ (x)Jν (y)|π 0 (~
p)i
is only suppressed as e−mπ Rmax . = Fµ,ν (x − y, p)h0|π 0 (y)|π 0 (~
p)i . (38)
We separately calculate the long-distance π 0 -exchange
contribution as illustrated in Fig. 8. We use the lattice The above relation suggest that the time-ordered prod-
calculation of the π 0 transition form factor: uct T Jµ (x)Jν (y), when used in between the vacuum
state and a single π 0 state, can be viewed as a properly
Fµ,ν (x, p) = h0|T Jµ (x)Jν (0)|π 0 (~
p)i (33) weighted π 0 interpolating field: Fµ,ν (x− y, p)π 0 (y). This
2 property can be used to calculate the following three-
= ǫµ,ν,ρ,σ xρ pσ F (x , p · x), (34)
point function:
where p is the Euclidean four-momentum of the π 0 state.
We use p = (imπ , ~0) = imπ t̂ in the numerical lattice hT Jµ (x)Jν (y)π 0 (z)i (39)
calculation. So we can obtain the numerical values of

For very large mπ |y − z| and y − z along the positive time
Fµ,ν x, imπ t̂ (35) direction, i.e. y − z = (|y − z|, ~0), the intermediate states
for all possible Euclidean space-time locations x directly between T Jµ (x)Jν (y) and π 0 (z) will be mostly π 0 states
from lattice QCD calculations, without assuming any with ~p ≈ ~0. We can employ Eq. (38) with p ≈ (imπ , ~0)
particular function form for the transition form factor. and obtain:
This form factor can then be used to construct approxi-
mations for infinite volume hadronic correlation functions hT Jµ (x)Jν (y)π 0 (z)i
 y−z 
with large separation. A detailed derivation is given in
≈ Fµ,ν x − y, imπ hT π 0 (y)π 0 (z)i. (40)
Appendix B. Here, we briefly describe the idea. First, we |y − z|
introduce a properly normalized π 0 interpolating field,
e.g. Note that we have expressed the above relation in an
O(4) rotationally covariant form, so the equation is also
i

π 0 (x) = Zπ1/2 √ ¯
ū(x)γ5 u(x) − d(x)γ5 d(x) . (36) valid for y − z in directions other than the time direc-
2 tion. Then, we can introduce the infinite-volume free
 11

48I strange con summand

48I strange discon summand
48I strange summand y
1.5

1 y′
x′
aµ × 1010 summand

0.5
x
0

−0.5

FIG. 8. The long-distance HLbL contribution to the muon

−1 g − 2 associated with π 0 exchange. The amplitudes inside the
boxes are calculated in lattice QCD while the pion propagator
linking them (solid lines) is given by the analytic, infinite-
−1.5
0 0.5 1 1.5 2 2.5 3 volume, continuum expression.
Rmax (fm)
48I strange con
48I strange discon
48I strange In the above approximation, |x − x′ | and |y − y ′ | will not
0.4
be too large in order to create localized pion states, so the
form factor F can be computed using a reasonably sized
0.2 lattice. As described in Eq. (33) and Eq. (35), we only
directly calculate the transition form factor F for a zero-
aµ × 1010 partial-sum

0 momentum pion state. To obtain the form factor needed

in Eq. (42) above, we need to perform O(4) rotations:
−0.2    
Fµ,ν x e, imπ n̂ = Λµ,µ′ Λν,ν ′ Fµ′ ,ν ′ xe′ , imπ t̂ , (43)
−0.4
where x̃ can be either x′ −x or y ′ −y in Eq. (42) and n̂ can
be ±(x − y)/|x − y|. The O(4) rotation matrix Λ = Λ(n̂)
−0.6
and Euclidean space-time coordinate x e′ satisfy:

−0.8 δµ,ν = Λµ,µ′ Λν,ν ′ δµ′ ,ν ′ , (44)

0 0.5 1 1.5 2 2.5 3
Rmax (fm) n̂µ = Λµ,µ′ t̂µ′ , (45)
eµ = Λµ,µ′ x
x e′µ′ . (46)
FIG. 7. Strange quark contributions computed on the 48I
ensemble from the connected diagrams, the disconnected dia- We would like to emphasize again that we do not as-
grams, and the total. The upper plot shows the corresponding sume any particular form of the π 0 transition form fac-
summands and the lower plot shows the partial sum. tor and perform fits. The inputs to Eq. (42) can be di-
rectly obtained from the O(4) rotation in Eq. (43) and
the lattice QCD calculation of the position-space matrix
scalar propagator Dπ0 (z), with physical pion mass. For elements Fµ,ν (x, imπ t̂) for a zero-momentum π 0 state de-
sufficiently large |x − y|, we have: fined in Eq. (33) and Eq. (35). The only approximation
hT π 0 (x)π 0 (y)i = Dπ0 (x − y). (41) is made in Eq. (42), which is based on large mπ |x − y|.
Using the 48I ensemble to calculate F and making this
Now, the relation between Eq. (33) and Eq. (40) is estab- approximation, we obtained the long-distance part cor-
lished. We can approximate the infinite volume hadronic responding to Rmax > 4 fm:
four-point function illustrated in Fig. 8 in the large |x−y|
region in a similar approach: aµ (Rmax > 4 fm) × 1010 (47)
0
hT Jµ′ (x′ )Jµ (x)Jν ′ (y ′ )Jν (y)i ≈ Dπ0 (x − y) ≈ aµπ -exch
(Rmax > 4 fm) × 1010 = 2.00(11)stat(28)syst ,
 x−y 
× Fµ′ ,µ x′ − x, imπ where we use the same subtracted infinite-volume QED
|x − y|
 weighting function as in the direct calculation. We es-
y−x  timate the relative systematic uncertainty from the ap-
× Fν ′ ,ν y ′ − y, imπ . (42)
|y − x| proximation employed in Eq. (42) due to not including
 12

the charged loop contribution as: since the on-shell condition p2 = −m2π for the physical
π 0 state cannot be satisfied always in the inverse Fourier
e−2mπ (4 fm) transformation. With the form factor defined above, we
≈ 6% (48)
e−mπ (4 fm) can construct the π 0 -pole contribution to the hadronic
four-point function:
and the error from approximating the π 0 propagation
direction to be along x − y to be: hT Jµ (x)Jµ′ (x′ )Jν (y)Jν ′ (y ′ )i (53)
Z Z
0.5 fm ≈ d u d4 v Dπ0 (u − v)
4
≈ 13% (49)
4 fm 
× Feµ,µ′ (u, x, x′ )Feν,ν ′ (v, y, y ′ )
where we assume the typical separation for |x − x′ | and
|y − y ′ | to be 0.5 fm and |x − y| to be 4 fm. Combin- + Feµ,ν (u, x, y)Feµ′ ,ν ′ (v, x′ , y ′ )
ing these two estimates in quadrature, we assigned 14% 
total systematic uncertainty for the long distance part + Feµ,ν ′ (u, x, y ′ )Feµ′ ,ν (v, x′ , y)
(Rmax > 4 fm) contribution. Using the ratio between the
connected and disconnected contributions as described in where Dπ0 (x − y) is the infinite volume free scalar propa-
Eq. 28, we can also properly assign this contribution to gator with physical pion mass. In this calculation, we use
the connected and disconnected diagrams. The results the Lowest Meson Dominance (LMD) model [38, 67, 68]
are shown in Tab. II with the labels “. . . Rmax > 4fm”. for the π 0 transition form factors.
Fπ0 γγ (q12 , q22 ) ≈ FπLMD 2 2
0 γγ (q1 , q2 )

D. Finite volume corrections 8π 2 Fπ2

=FπVMD 2 2
0 γγ (q1 , q2 ) + (54)
3m2V
 
The long-distance part of the HLbL contribution usu-
× FπTE 2 2 VMD 2 2
0 γγ (q1 , q2 ) − Fπ 0 γγ (q1 , q2 ) ,
ally suffers the most significant finite volume effects.
However, the long-distance π 0 exchange contribution cal-
where
culated in Section IV C is performed in infinite volume
and is free of finite volume effects. Therefore, we only m2V m2V
FπVMD 2 2
0 γγ (q1 , q2 ) = (55)
need to correct for the finite volume effects for the rela- q12 + mV q2 + m2V
2 2
tively short-distance region, where Rmax < 4 fm. These
finite volume effects are expected to be quite small due
to the constraint Rmax < 4 fm. The finite volume effects m2V /2 m2V /2
FπTE 2 2
0 γγ (q1 , q2 ) = + . (56)
can be estimated with the “π 0 -pole” contribution as de- q12 + m2V q22 + m2V
fined in Ref. [20]. Similar to Ref. [21, 66], we define the
The parameters used in the calculation are mπ =
π 0 transition form factor in Euclidean space-time as:
134.9766 MeV, Fπ = 92 MeV, mV = 770 MeV. We
Fµ,ν (x, p) = h0|T Jµ (x)Jν (0)|π 0 (~ p)i discretize the model and calculate it with lattices of dif-
Z 4 ferent sizes using the same subtracted infinite volume
d q1 iq1 ·x −i QED weighting as the direct calculation. The results are
= e ǫµ,ν,ρ,σ q1 ρ q2 σ Fπ0 γγ (q12 , q22 ) (50)
(2π)4 4π 2 Fπ given in Tab. V. Note the physical size of the “Match”
and “Coarse” lattices are the same as the 48I ensemble,
where p = q1 + q2 and p2 = −m2π . We can convert the which we used to perform our main lattice QCD calcula-
form factor into coordinate space via the following: tion. The spatial size is L = 5.5 fm. The “Large” lattice
Z has the same lattice spacing as the “Coarse” lattice but
d4 u eip·u Feµ,ν (u, x, y) = Fµ,ν (x − y, p). (51) has a much larger physical size, L = 16.4 fm. We use the
difference between the “Large” and the “Coarse” results
of aµ (Rmax < 4 fm)×1010 as the finite volume correction.
The solution to the above condition is not unique. Based
Comparing the results of “Match” and “Coarse”, we esti-
on the momentum space form factor Fπ0 γγ (q12 , q22 ) in
mate the uncertainty due to the non-zero lattice spacing
Eq. (50), we obtain
effects of this model lattice calculation can be about 12%.
i Combined with the additional 20% uncertainty due to the
Feµ,ν (u, x, y) = ǫµ,ν,ρ,σ ∂ρx ∂σy (52) inaccuracy of the model itself (and potential finite vol-
4π 2 Fπ
Z 4 Z 4 ume corrections of other heavier intermediate states), we
d q1 iq1 ·(x−u) d q2 iq2 ·(y−u) obtain our final estimate of the finite volume correction:
× e e Fπ0 γγ (q12 , q22 ).
(2π)4 (2π)4
aFV-corr
µ (Rmax < 4 fm) × 1010 = −0.47(11)syst. (57)
Note this definition of the position space form factor
is different from the inverse Fourier transformation of Similar to the long-distance π 0 -exchange contribution, we
the Euclidean space hadronic matrix element Fµ,ν (x, p), use the ratio between the connected and disconnected
 13

contributions as described in Eq. 28 to assign this contri- correction with this form instead, the size of the correc-
bution to the connected and disconnected diagrams. The tion is 0.20(4) × 10−10 . In Ref. [69], the HLbL is calcu-
results are shown in Tab. II as “48I light con FV-corr”, lated with Chiral effective theory, we have calculated the
“48I light discon FV-corr”, and “48I light FV-corr”. pion mass dependence with Eq. (13) of this reference.
Also, note that we can calculate the long-distance π 0 - The shift of aµ due to the pion mass change from the
pole contribution using this LMD model. In Tab. V, 139 MeV used in the lattice calculation to the physical
we also list the contribution to the region Rmax > 4 fm. point 135 MeV is 0.25 × 10−10 . Here, we do not include
Using the LMD model results calculated with the “Large” the correction to the charged pion loop contribution of
lattice, we obtain aµ due to the pion mass shift, as the charged pion loop
0
contribution are significantly reduced compared with the
aµπ -pole;LMD
(Rmax > 4 fm) × 1010 = 2.34. (58) scalar QED estimate due to the pion dressing [70–72].
We, therefore, estimate a 50% systematic uncertainty to
This result agrees with the long-distance π 0 -exchange these corrections to account for possible inaccuracies of
contribution calculated with π 0 transition form factors the empirical chiral extrapolations, plus any systematic
from the previously described lattice QCD calculation effects in obtaining the 32D and 24DH results, including
and long-distance approximation. the effects caused by the hybrid method used to calcu-
late the long-distance part of the disconnected diagrams,
the discretization effects, the finite volume effects, and so
E. Pion mass extrapolation on. The final results for the correction due to the slight
mismatch of the pion mass is:
While the calculation is performed very near the phys-
ical pion mass, there is a slight mis-tuning of the lat- aµmπ -corr × 1010 (61)
tice parameters which leads to a slightly heavier mπ =
= (aµ (mπ = 135 MeV) − aµ (mπ = 139 MeV)) × 1010
139 MeV compared to the physical pion mass mπ =
134.9766 MeV. We correct this small difference using the = 0.35(0.07)stat(0.17)syst .
32D and 24DH ensembles. These two ensembles have the
same lattice spacing but different pion masses. The re- We also show the results in Tab. II as “48I light con
sults of these two ensembles are listed in Tab. VI. Note mπ -corr”, “48I light discon mπ -corr”, and “48I light mπ -
that for the long-distance contribution (Rmax > 4 fm) for corr”.
the 32D physical pion mass ensemble, we use the same
long-distance π 0 -exchange contribution calculated with
F. Continuum limit extrapolation
ensemble 48I as described in Section IV C. Similar to the
48I calculation, we calculate the “no-pion” contribution
and obtain the “hybrid-2fm” results for the “discon” and In our previous work [24], we used finite-volume lat-
“total” results. To reduce the statistical error, we use tices for the entire calculation, including QED. The
the “hybrid-2fm” results to calculate the pion mass cor- QEDL scheme was used for the photon propagators, and
rection. To obtain the pion mass correction to our main sizable discretization effects were observed. We took the
48I results, we assume the following pion mass depen- continuum limit with the 48I and 64I ensembles. These
dence for the “no-pion” contribution: two ensembles have different lattice spacings but all other
properties are almost identical. They have the same
ano-pion
µ (mπ ) = C1 + C2 m2π . (59) gauge and fermion actions, and almost the same pion
mass and volume in physical units. We used the results
Due to the pion pole piece in the quark-connected dia- from these two ensembles to perform the continuum ex-
gram, we assume a more singular pion mass dependence trapolation. The continuum limit value is 30(12)% higher
for the “con” contribution: compared to the value on the 48I ensemble for the con-
C2 nected piece and 58(27)% higher for the disconnected
acon
µ (mπ ) = C1 + . (60) piece. The numbers in parentheses are statistical un-
m2π
certainties.
The disconnected diagram contribution can be obtained In this work, we have mainly used the 48I ensemble
as a combination of the “con” and “no-pion” contribution to perform the calculation with the infinite volume, con-
as given in Eq. (31)((32)). tinuum QED weighting function described in Ref. [46].
The fitting form in Eq. (60)(plus an additional m2π Compared with QEDL , the infinite volume QED weight-
term) was used for the chiral extrapolation of the con- ing function reduces the finite volume effects from power-
nected diagram in Ref. [48] and seems to provide the most law to exponential.
accurate description of the connected diagram’s mass One may expect similar discretization effects for these
dependence. The correction from the chiral extrapola- two calculations since they should be the same in the
tion is 0.35(7) × 10−10 with the above fitting form. In large volume limit and only differ by finite volume ef-
Ref. [26], the π 0 -pole contribution is calculated to be- fects. However, there are two reasons that the discretiza-
have as C1 + C2 log2 (m2π ), If we estimate the pion mass tion errors in the current work should be much smaller.
 14

Label Size a−1 /GeV aµ (Rmax < 4 fm) × 1010 aµ (Rmax > 4 fm) × 1010
Match 483 × 96 1.73 5.19 0.22
Coarse 243 × 48 0.865 4.65 0.20
Large 643 × 128 0.865 4.18 2.34

TABLE V. Lattice calculation of the π 0 -pole contribution using the LMD model. Note that this is not a lattice QCD calculation.

name Rmax limit aµ (32D) ×1010 aµ (24DH) ×1010

con 2 fm 8.17 (0.32)stat 8.29 (0.10)stat
con 4 fm 21.29 (3.37)stat 12.35 (0.24)stat
con ∞ 28.84 (3.39)stat 12.35 (0.24)stat
no-pion 2 fm 3.91 (0.63)stat 4.18 (0.57)stat
no-pion 4 fm 4.96 (4.46)stat 6.74 (4.18)stat
discon 2 fm −2.10 (0.59)stat −1.91 (0.55)stat
discon 4 fm −10.70 (3.84)stat −2.34 (4.18)stat
discon ∞ −16.25 (3.84)stat −2.34 (4.18)stat
discon hybrid-2fm 4 fm −11.74 (2.52)stat −4.90 (0.54)stat
discon hybrid-2fm ∞ −17.30 (2.53)stat −4.90 (0.54)stat
total 2 fm 6.07 (0.66)stat 6.38 (0.57)stat
total 4 fm 10.59 (4.97)stat 10.01 (4.18)stat
total ∞ 12.59 (4.98)stat 10.01 (4.18)stat
total hybrid-2fm 4 fm 9.55 (1.12)stat 7.45 (0.59)stat
total hybrid-2fm ∞ 11.55 (1.12)stat 7.45 (0.59)stat

TABLE VI. Study of the pion mass dependence using the 32D and 24DH ensembles. For the 32D (physical pion mass
ensemble) results, we apply the 48I long-distance π 0 -exchange contribution in the Rmax > 4 fm region. For the 24DH ensemble
(mπ = 341 MeV), we assume the contribution is negligible in the Rmax > 4 fm region.

First, the infinite volume QED weighting function is cal- The study of the continuum limit of the strange quark
culated in the continuum with a semi-analytic method, connected diagrams strongly suggests the light quark
which eliminates the discretization errors in the QEDL connected diagrams, in particular in the relatively short
weighting function which is calculated on a discrete lat- distance region where we performed a reliable continuum
tice. Second, and more importantly, we apply the sub- extrapolation for the strange quark connected diagrams,
traction scheme as described in Ref. [46] and repeated is under 8%. We also have the same 64I dataset for the
in Eq. 12. In a test calculation of the leptonic light-by- light quark connected diagram. We have listed the com-
light contribution to the muon g − 2, we found about a parison in Tab. VII. The results from the 64I ensemble
factor of 5 reduction in the discretization error after this are statistically consistent with the 48I ensemble. Un-
subtraction is applied. We expect a similar reduction forunately, the statistical error for 64I is relatively larger
in discretization error in the hadronic case, and this is due to the absence of the improvement we made for the
indeed the case in our calculation of the strange quark new 48I calculation. We cannot constrain the continuum
connected diagrams. extrapolation for the light quark connected diagram bet-
The strange quark connected contribution is shown in ter than the 8% estimate above based on the strange
Fig. 9. Compared with the light quark, due to the heav- quark connected results.
ier strange quark mass, the statistical noise in the long-
In the long distance (large Rmax ) region, the contribu-
distance region is very small. The overall signal-to-noise
tions are mostly coming from the π 0 -exchange contribu-
ratio is much better than the light quark case which al-
tion. As can be seen from Eq. (42), the hadronic four-
lows a precise continuum limit extrapolation. As is shown
point function in the HLbL diagrams in the long-distance
in the figure, the limit is 7.9(2.6)% higher than the 48I
part can be approximated by the product of two such
results. The continuum limit is obtained by performing
transition form factors. Therefore, the discretization ef-
an O(a2 ) extrapolation to a → 0,
fects of the direct calculation of the four-point function
aµ (a2 ) = C1 + C2 a2 . (62) in the long distance region will be similar to the dis-
cretization effects of the calculation of the square of the
Compared with the observed ∼ 30% discretization error transition form factor. In the π 0 → e+ e− calculation of
in the QEDL calculation of the light quark connected Ref. [73], this relevant π 0 → γγ transition form factor is
diagrams, there is indeed a significant reduction in the studied with the same 48I ensemble as this calculation.
relative discretization error with the subtracted infinite In addition, the continuum limit is calculated with a par-
volume weighting function, despite the heavier strange allel 64I ensemble calculation. We quote the corrections
quark mass, which usually leads to larger discretization obtained in the calculation in Tab. VIII. Note that the
errors. correction is zero consistent given the statistical uncer-
 15

Observable 48I 64I Percentage correction

acon
µ (Rmax < 1.0 fm) × 1010
1.27(0.02) 1.32(0.05) 9.2(8.6)
acon
µ (Rmax < 1.5 fm) × 10
10
3.92(0.05) 4.04(0.23) 7.0(12.6)
aµ (Rmax < 2.0 fm) × 1010 7.24(0.15)
con
7.19(0.66) −1.4(20.0)
acon
µ (Rmax < 2.5 fm) × 10
10
10.64(0.31) 9.61(1.40) −20.8(29.2)

TABLE VII. The last column is the percent correction to obtain the continuum limit relative to the value of the 48I ensemble.
The values in parentheses are the statistical uncertainty of the corresponding quantity.

48I strange con summand Observable Percentage correction

64I strange con summand (ImA)2 -11.4(11.6)%
strange con summand
0.35 (ReA)2 -16.8(15.8)%

0.3 TABLE VIII. Percent correction to obtain the continuum

limit relative to the value computed on the 48I ensemble in the
0.25 π 0 → e+ e− calculation in Ref. [73]. The values in parentheses
are the statistical uncertainty of the corresponding quantity.
aµ × 1010 summand

0.2 The square of the amplitude is used to estimate the lattice

spacing error of the long-distance π 0 -exchange contribution
0.15 to the muon g − 2 in this work.

0.1
For the light quark disconnected diagrams, we again
0.05 separately discuss the short and long distance regions.
In Fig. 5, we observe that in the short distance region
0 (Rmax . 2 fm), the disconnected diagram contribution is
much smaller than the connected contribution. This is
−0.05
0 1 2 3 4 5 expected as the disconnected diagrams involve two quark
Rmax (fm) loops and the exchange of at least two gluons. Therefore,
48I strange con the disconnected diagrams are suppressed by O(α2s ) and
64I strange con 1/Nc relative to the connected diagrams. Due to the
strange con smallness of the contributions from the short-distance
0.4
disconnected diagrams, we expect that the same reason-
0.35 ing applies to the discretization effects.
For the long distance region, we again have Eq. (28),
0.3
which mandates the relationship with the connected dia-
aµ × 1010 partial-sum

0.25 grams, even for non-zero lattice spacing and finite volume
(so long that the lattice is large enough to contain the
0.2
correlation function). Therefore, we expect the long dis-
0.15 tance part of the disconnected diagrams to have the same
relative discretization effects as the connected diagrams
0.1
and also the combined total. Note that the situation is
0.05 different from the Mainz group’s work in Ref. [47, 48].
While the argument presented here should also apply
0 to the hadronic function in their calculation, the differ-
−0.05 ent QED kernels (due to the subtraction scheme of the
0 1 2 3 4 5 QED kernel and the choice of xref ) used for the connected
Rmax (fm) and disconnected diagrams lead to different shapes of the
summands and also different discretization errors.
FIG. 9. Strange quark connected part of aHLbL
µ (lower panel). Without a 64I ensemble calculation of similar preci-
Continuum limit (upper curve) and finite lattice spacing re- sion, the continuum limit cannot be reliably taken. How-
sults on 48I and 64I ensembles. Corresponding summands ever, based on the evidence presented above, in particular
(upper panel). the size of the correction in the strange quark connected
diagrams, we estimate an 8% uncertainty from non-zero
lattice spacing for the contributions from the light quark
connected, disconnected, and total contributions, and the
tainty. In Fig. 4 of Ref. [73], the partial sum of these strange quark disconnected diagrams. The corrections
quantities for both 48I and 64I are plotted. More sta- for light quark are shown in Tab. II as “48I light con
tistically precise agreement is observed at smaller time a2 -corr”, “48I light discon a2 -corr”, and “48I light a2 -
separation of the two electromagnetic vector currents. corr”. The corrections for the strange quark are shown in
 16

Tab. IV as “48I strange con a2 -corr”, “48I strange discon contribution is included as 0.3(0.1) × 10−10 , calculated
a2 -corr”, and “48I strange a2 -corr”. based on perturbation theory and consideration of charm
meson resonances [23]. The value is similar to the strange
quark connected diagram contribution due to the heavier
G. Sub-leading disconnected diagrams mass of the charm quark and larger electric charge. In
a recent work by the Mainz group [49], the charm quark
So far, our discussions have centered on the first two contribution is calculated with lattice QCD, including
diagrams in Fig. 2. The remaining diagrams, which are both the charm quark connected and disconnected dia-
expected to be small due to the suppression by flavor grams. The result for the total charm quark contribution
SU(3) and quark electric charge factors, are not yet in- is 0.28(0.05)×10−10, where the uncertainty is mostly due
cluded. In our earlier work [24], we have already calcu- to the systematic effects from modeling the lattice spac-
lated the third diagram in Fig. 2, which is expected to ing and Mηc dependence. In this work, we use this more
be the largest within the remaining the sub-leading dis- recent lattice calculated result to account for the contri-
connected diagrams based on the flavor SU(3) and quark bution from the charm quark. The values are shown in
electric charge factor counting. To estimate the contribu- Tab. IX as “charm con”, “charm discon” “charm”.
tion from this diagram, we employ the same subtracted
infinite volume QED weighting function and the 24D en-
semble for the quark loops. The result, shown in Fig. 10, V. CONCLUSION
is zero consistent, and we estimate this contribution to
be 0.0(0.5) × 10−10 . We summarized all the results discussed above in
Tabs. II, IV, and IX. Adding all the individual contri-
24D sub-leading-discon butions and corrections, we obtain our final result for
0.4 the HLbL contribution to the muon g − 2:

0.2 aHLbL
µ × 1010 = 12.47(1.15)stat(0.95)syst [1.49], (63)

0 where the systematic errors were added in quadrature.

aµ × 1010 partial-sum

We also give the HLbL contributions to the muon g − 2

−0.2 from the connected and disconnected pieces separately:

−0.4
aHLbL,con
µ × 1010 = 26.36(1.33)stat(1.99)syst [2.39],
(64)
−0.6 aHLbL,discon
µ × 1010 = −13.89(1.47)stat(1.22)syst [1.91].
(65)
−0.8
Numbers in square square brackets denote total erorrs,
−1 combining the statistical uncertainty (“stat”) and sys-
0 0.5 1 1.5 2 2.5 3 tematic ones (“sys”) in quadrature.
Rmax (fm) Here, we summarize techniques used in this calculation
that we believe are important to obtain the precision for
FIG. 10. Sub-leading disconnected diagrams computed on the HLbL scattering at the physical pion mass:
24D ensemble. [24]
• The subtracted infinite volume QED weighting
In the more recent work by the Mainz group [48], all function developed in our previous work [46].
the sub-leading diagrams have been calculated. The re-
sult is still zero consistent but a more stringent bound • Use all combinations of the point source propaga-
is obtained. Therefore, in this work, we use their result tors to calculate the HLbL diagrams as described
to account for the contribution from all the sub-leading in section III.
disconnected diagrams. The value is shown in Tab. IX
• The adaptive sampling scheme used for the dis-
as “sub-leading discon”.
connected diagram calculation as described in sec-
tion III.
H. Charm quark contributions • The rearrangement of the connected and discon-
nected diagrams in Eq. (29,32) based on Eq. (28).
The charm quark contribution is also expected to be
very small due its heavy mass relative to the light quark • The AMA method [58] and efficient GPU solver for
and the strange quark. In the 2020 white paper [5], its MDWF propagators from Grid and GPT.
 17

Contribution aµ ×1010
light con 25.70 (1.33)stat (1.99)syst [2.39]
light discon −12.71 (1.87)stat (1.17)syst [2.20]
light discon hybrid-2.5fm −13.50 (1.36)stat (1.21)syst [1.82]
light discon hybrid-2fm −13.37 (1.29)stat (1.46)syst [1.95]
light total 12.99 (2.11)stat (0.90)syst [2.29]
light total hybrid-2.5fm 12.20 (1.01)stat (0.95)syst [1.38]
light total hybrid-2fm 12.33 (0.90)stat (1.25)syst [1.55]
strange con 0.35 (0.01)stat
strange discon −0.38 (0.61)stat (0.03)syst [0.61]
strange discon hybrid-2.5fm −0.36 (0.22)stat (0.03)syst [0.22]
strange discon hybrid-2fm −0.28 (0.12)stat (0.04)syst [0.13]
strange total −0.03 (0.61)stat (0.03)syst [0.61]
strange total hybrid-2.5fm −0.00 (0.22)stat (0.03)syst [0.23]
strange total hybrid-2fm 0.07 (0.12)stat (0.04)syst [0.13]
sub-leading discon 0.00 (0.07)syst [48]
charm con 0.31 (0.04)syst [49]
charm discon −0.03 (0.02)syst [49]
charm total 0.28 (0.05)syst [49]
con 26.36 (1.33)stat (1.99)syst [2.39]
discon −13.12 (2.30)stat (1.18)syst [2.59]
discon hybrid-2.5fm −13.89 (1.47)stat (1.22)syst [1.91]
discon hybrid-2fm −13.68 (1.35)stat (1.47)syst [1.99]
total 13.24 (2.53)stat (0.90)syst [2.68]
total hybrid-2.5fm 12.47 (1.15)stat (0.95)syst [1.49]
total hybrid-2fm 12.68 (0.98)stat (1.26)syst [1.59]

TABLE IX. Summary of final results. Values are different contributions to aµ ×1010 where aµ = (gµ −2)/2. The numbers in the
square bracket (except references) are the statistical and systematic uncertainty combined in quadrature. The “strange discon”
contribution includes disconnected diagrams where one or two loops are strange quark loops. The tag “hybrid-” indicates the
same quantity obtained using the hybrid approach as described in Section IV A and IV B. The result for “total hybrid-2.5fm”
is used as our final result.

aµ ×1010 Hadron Models

all con QEDL 24.46 (2.35)stat (5.11)syst [5.62] Dispersive + Data
all con diff 1.90 (2.76)stat (5.48)syst [6.14] Lattice
all discon QEDL −16.45 (2.09)stat (3.99)syst [4.50]
all discon diff 2.56 (2.57)stat (4.17)syst [4.90] PdRV 09
total QEDL 8.17 (3.03)stat (1.77)syst [3.51] N/JN 09
total diff 4.30 (3.25)stat (2.01)syst [3.82]
FJ 17
TABLE X. Comparison of this work with our previous QEDL White paper 20
results [24] plus the charm quark contribution from [5, 23].
The field “all con” includes all connected diagram contribu- RBC-UKQCD 19
tions, including the charm quark. The field “all discon” in-
Mainz 21
cludes all the disconnected diagram contributions, including
the sub-leading disconnected diagrams and disconnected dia- This Work
grams including charm quarks. The fields “all con diff”, “all
discon diff”, “total diff” shows the difference obtained by sub- 4 6 8 10 12 14 16 18
tracting the QEDL results from the new results presented in
this work. The statistical and systematic uncertainties of the aHLbL
µ
10
× 10
two works are almost independent. We therefore add them
in quadrature for the difference. The numbers in the square FIG. 11. Comparison of our result with values in the litera-
bracket are the statistical and systematic uncertainty com- ture. The hadronic model values are from Refs. [26–29]. The
bined in quadrature. dispersive data driven result is compiled in Ref. [5]. The lat-
tice results include Refs. [24, 48, 49] and this work.

We can compare these results with our previous

work [24] based on the finite volume QEDL formulation. rent result is 1.12 standard deviations higher than the
The comparison is summarized in Table X. We can see previous results, possibly due to a slightly larger statis-
that the results for both the connected and disconnected tical fluctuation. We also compare the final result in this
diagrams are in good agreement. For the total, the cur- work with the existing literature in Fig. 11. The new
 18

result is consistent with previous determinations. The γµ matrices satisfy the Euclidean space-time metric
γµ γν + γν γµ = 2δµ,ν , (A5)
ACKNOWLEDGMENTS and
γ5 = γx γy γz γt . (A6)
We thank our RBC and UKQCD collaborators for
helpful discussions and critical software and hardware
support. T.B. and L.J. have been supported under
Appendix B: π 0 long-distance contribution
US DOE grant DE-SC0010339. L.J. is also supported
by US DOE Office of Science Early Career Award DE-
SC0021147. N.C. is supported by US DOE grant DE- In Eq. (42), we replaced the QCD, Euclidean space-
SC0011941. M.H. is supported by Japan Grants-in- time, four-current connected Green’s function with the
Aid for Scientific Research, No20K03926. T.I., C.J., product of two amplitudes, each coupling a pair of cur-
and C.L. were supported in part by US DOE Contract rents to an on-shell π 0 . These two amplitudes are joined
DESC0012704(BNL). T.I. and C.J. were supported in by a pion propagator and all amplitudes are expressed
part by the Scientific Discovery through Advanced Com- in position space, so they can be directly inserted in our
puting (SciDAC) program LAB 22-2580. T.I. is also standard position-space evaluation of the HLbL ampli-
supported by the Department of Energy, Laboratory tude. Since the final expression involves two indepen-
Directed Research and Development (LDRD No. 23- dent factors evaluating the πγγ coupling which are con-
051) of BNL and RIKEN BNL Research Center, and nected by an analytic, position-space pion propagator,
by JSPS KAKENHI under grant numbers JP26400261, this QCD part of the HLbL amplitude can be evaluated
JP17H02906. C.L. has been supported by a DOE Of- in a “QCD volume” of arbitrary size. In particular, this
fice of Science Early Career Award. We developed volume could be much larger than that of the gauge con-
the computational code used for this work based on figurations used to compute each πγγ vertex function.
the Columbia Physics System (CPS), Grid, GPT, and Here we work out a concrete derivation of this formula
Qlattice. Computations were performed mainly under that can be used to evaluate the long-distance part of
the ALCC Program of the US DOE on the SUMMIT the π 0 exchange contribution to the leading order in 1/L.
computer at the Oak Ridge Leadership Computing Fa- We leave open the possibility that this approach could be
cility, a DOE Office of Science Facility supported under developed further to systematically capture terms falling
Contract No. DE-AC05-00OR22725. These calculations with higher powers of 1/L if the large volume π 0 con-
used gauge configurations and propagators created using tribution is expressed as a power series in 1/Ln where
resources of the Argonne Leadership Computing Facility, L is the size of the QCD volume. We assume the QED
which is a DOE Office of Science User Facility supported volume to be infinite.
under Contract DE-AC02-06CH11357. We also acknowl- We begin with the Euclidean-space Green’s function
edge computer resources at the Oakforest-PACS super- defined in Eq. (15):
computer system at Tokyo University, partly through the 6e4 Hµ′ ,µ,ν ′ ,ν (x′ , x, y ′ , y) (B1)
HPCI System Research Project (hp180151, hp190137), ′ ′


the BNL SDCC computer clusters at Brookhaven Na- = T Jµ′ (x )Jµ (x)Jν ′ (y )Jν (y) .
tional Lab as well as computing resources provided We will choose x′ and x close to each other as are y ′ and
through USQCD at Brookhaven and Jefferson National y. However, we will assume that the (x′ , x) pair is far
Labs. from the (y ′ , y) pair. Define:
e = x′ − x,
x (B2)
Appendix A: Notation ye = y ′ − y. (B3)
Since Euclidean space-time is rotational invariant, with-
We use Sµ and G to denote free muon and photon
out loss of generality, we can choose x − y to be along the
propagators:
Euclidean time direction, with x − y = (t, ~0) and t > 0.
Z We now follow the usual steps to obtain a variant of
d4 p 1
Sµ (x, y) = eip·(x−y) (A1) the Källen-Lehman representation of the amplitude but
(2π)4 i 6 p + m
Z keep only the single π 0 intermediate state since, as the
d4 p 1 lightest particle, its exchange will dominate this Green’s
= (− 6 ∂ x + m) eip·(x−y) , (A2)
(2π)4 p2 + m2 function when x and y are far separated:
0
Z 6e4 Hµπ′ ,µ,ν ′ ,ν (x′ , x, y ′ , y) (B4)
d4 p 1 ip·(x−y) Z
G(x, y) = e (A3) 3
d p 1

(2π)4 p2 = 0 T Jµ′ (x′ )Jµ (x) π 0 (~ p)
1 1 (2π)3 2Eπ,~p
= . (A4)

4π 2 (x − y)2 × π 0 (~ p) T Jν ′ (y ′ )Jν (y) 0 ,
 19

where the superscript π 0 indicates that only the π 0 con- order in 1/|x − y|:
tribution to H is represented.
YN  
Next we use four-dimensional translational invariance ∂
Dπ0 (x − y)
to remove the variables x and y from the two current- i=1
∂xρi
current-π 0 amplitudes. We will also replace these two N   −mπ |x−y|
O(4)-covariant current-current-π 0 amplitudes by func- Y ∂ ce
= (B10)
tions of the four momentum p = (iEπ,~p , ~ p). Given the
i=1
∂x ρi |x − y|3/2
factors of exp (ip · x) and exp (−ip · y) which appear, we N  
can then replace the four vector pµ by −i∂/∂xµ or i∂/∂yµ Y (x − y)ρi c e−mπ |x−y|
≈ −mπ , (B11)
as needed. For example, |x − y| |x − y|3/2
i=1

where the relation holds for arbitrary N and ρi labels

0 T Jµ′ (x′ )Jµ (x) π 0 (~
p)

the component of x appearing in the ith derivative in the
= 0 T Jµ′ (e x)Jµ (0) π 0 (~
p) eip·x product of N derivatives. Here we use the large distance
x, p)eip·x
= Fµ′ ,µ (e behavior of Dπ0 (x − y):
 
∂ c e−mπ |x−y|
= Fµ′ ,µ x e, −i eip·x . (B5) Dπ0 (x − y) ≈ . (B12)
∂x |x − y|3/2

where the definition of the π 0 transition form factors fol- Introducing this approximation into Eqs. (B9), we obtain
lows Eq. (33). 0
6e4 Hµπ′ ,µ,ν ′ ,ν (x′ , x, y ′ , y) (B13)
Using this approach to remove the explicit dependence
of the two amplitudes on pµ , we can then perform the ≈ Fµ′ ,µ (e
x, imπ n̂) Fν ′ ,ν (e
y , −imπ n̂) Dπ0 (x − y).
final step of the usual Källen-Lehman derivation and in-
where n̂µ is the Euclidean unit vector (x − y)µ /|x − y|.
troduce the free pion propagator:
This is finally the equation of interest, Eq. (42), which
Z appeared earlier. Since the two quantities expressed in
d3 p 1 terms of the form factor F can be evaluated from inde-
3
eip·(x−y) pendent lattice QCD ensemble averages, the expression
(2π) 2Eπ,~p
Z on the right hand side of Eq. (B13) can be evaluated from
d3 p 1
= ei~p·(~x−~y) e−Eπ,~p(xt −yt ) (B6) a standard lattice QCD ensemble average followed by an
(2π)3 2Eπ,~p O(4) rotation as defined in Eq. (43), which we repeat
Z
d4 p eip·(x−y) here:
= (B7)    
(2π)4 p2 + m2π
Fµ′ ,µ xe, imπ n̂ = Λµ′ ,ρ′ Λµ,ρ Fρ′ ,ρ xe′ , imπ t̂ , (B14)
= Dπ0 (x − y) (B8)
The O(4) rotation matrix Λ = Λ(n̂) and Euclidean space-
where Dπ0 (x − y) is the free Euclidean-space propagator time coordinate x̃′ satisfy:
for a scalar particle of mass mπ .
δµ,ν = Λµ,µ′ Λν,ν ′ δµ′ ,ν ′ , (B15)
We can then rewrite Eq. (B4) in terms of F and Dπ0
to obtain the result: n̂µ = Λµ,µ′ t̂µ′ , (B16)
′
eµ = Λµ,µ′ x
x eµ′ . (B17)
0
6e4 Hµπ′ ,µ,ν ′ ,ν (x′ , x, y ′ , y) (B9)  
    The rotated form factor Fρ′ ,ρ x e′ , imπ t̂ follows the def-
∂ ∂
= Fµ′ ,µ x e, −i Fν ′ ,ν ye, −i Dπ0 (x − y). inition Eq. (B5) or Eq. (33) with zero momentum π 0
∂xµ ∂yµ
state:
 
At this step, we have made explicit the O(4)-covariance Fρ′ ,ρ xe′ , imπ t̂
of the right-hand side of the equation. So this equation
should hold for an arbitrary orientation of x − y. = h0|T Jρ′ (e p = ~0)i
x′ )Jρ (0)|π 0 (~
So far all of the steps taken have been exact. We expect p x′ )Jρ (0)π 0 (−tsep )i
hT Jρ′ (e
= 2mπ L3 lim p , (B18)
the quantity computed, the contribution of a single pion ttsep →+∞ hT π 0 (tsep )π 0 (−tsep )i
exchange, to dominate the long distance limit in which
|x − y| is large, i.e. |x − y| & L. Now we will evaluate where π 0 (t) is a zero momentum π 0 interpolating op-
the derivatives with respect to x and y which appear erator. We use Coulomb fixed gauge wall source pion
in these equations but keep only the leading term in an operator in our numerical lattice calculation. The other
expansion in powers of 1/L. For example, to leading factor Fν ′ ,ν (e
y , −imπ n̂) can be calculated similarly.
 20

The above procedure allows a direct calculation of the • fit7:

π 0 contribution to the long distance part the HLbL am- 6
Rmax
plitude without relying on any parameterization or mod- f (Rmax ) = A/fm 3
e−2mπ Rmax /(fm·GeV)
elling of the π 0 transition form factor. (Rmax + (C fm)3 )2
In this calculation, we employ this long-distance π 0 - (C7)
exchange approximation in the following region:
• fit8:
Rmax = max(|x − y|, |y ′ − y|, |x − y ′ |) ≥ 4 fm, (B19) 6
Rmax
f (Rmax ) = A/fm 3
e−mπ Rmax /(fm·GeV)
as defined in Eq. (27). This long-distance hadronic func- (Rmax + (C fm)3 )2
tion is then combined with infinite volume QED weight- (C8)
ing function for HLbL as illustrated in Fig. 8. Note that
we need to sum over all permutations of x, y, y ′ . where we set mπ = 139 MeV and mη = 550 MeV. For
“fit1”, “fit2”, “fit3”, we always constrain the B parameter
to be larger or equal to mπ .
Appendix C: Summand fitting form In Figs. 12,13,14,15,16,17, and 18, we show the result-
ing fits along with the data.
In section IV A, we have combined the disconnected di- We also apply these fit functions to different contri-
agrams and connected diagrams with a certain fraction butions, including the strange quark connected diagram,
multiplied to form ano-pion according to Eq. (29). We the light quark connected diagrams only, and the light
µ
have fitted the long distance behavior of this combina- quark disconnected diagrams only. Table XI summarizes
tion with Eq. (30) and the results are shown in Fig. 6. the fits. The strange quark connected fits are displayed
In this appendix, we will explore other possible fitting in Figs. 19,20,21 and, for the light quark connected dia-
forms and will also apply these fitting functions to other grams, in Figs. 22,23,24,25. And finally, the light quark
contributions. disconnected results are shown in Figs. 26,27,28,29. The
The fitting forms are first fit function (fit1) defined in Eq. (C1) fits the contri-
bution from the strange quark connected diagrams, which
• fit1: This is the form used in the main calculation, are statistically very precise, very well. While the χ2d.o.f.
same as Eq. (29). for this fit is quite large, this is mostly due to the small
statistical error for the strange quark connected diagrams
6
Rmax results. The difference between the fit and data for these
f (Rmax ) = A/fm4 3
e−BRmax /(fm·GeV)
Rmax + (C fm)3 strange quark connected diagrams results is much smaller
(C1) than the needed precision for the fit for the light no-pion
combinations. Based on the comparison of qualities of
• fit2: all the fits performed, we chose fit1 defined in Eq. (C1)
as the central value for the tail part of the light quark
f (Rmax ) = A/fm4 Rmax
3
e−BRmax /(fm·GeV) (C2) no-pion contribution in our main calculation described in
section IV A.
• fit3: For the light quark no-pion long-distance contribu-
tions, we observe roughly 100% variation about the cen-
6
Rmax tral value (“fit1”), which is the reason for the 100% sys-
f (Rmax ) = A/fm 3
e−BRmax /(fm·GeV)
(Rmax + (C fm)3 )2 tematic uncertainty assigned to the contributions “48I
(C3) light no-pion Rmax > 2.5fm” and “48I light no-pion
Rmax > 2fm” in Table II. Note that the “fit3” result
• fit4: for “48I light no-pion” has large statistical uncertainty.
This is due to the fact that the data does not constrain
6
Rmax the exponent B very much, in contrast to “fit6” and
f (Rmax ) = A/fm4 3
e−mη Rmax /(fm·GeV)
Rmax + (C fm)3 “fit7” which describe the data well. Since the intermedi-
(C4) ate states should have at least energy 2mπ , we view the
results from “fit7” to be the upper bound for the fitting
• fit5: form “fit3”, which falls into the range of the results ob-
tained with “fit1” and the 100% systematic uncertainty
f (Rmax ) = A/fm4 Rmax
3
e−mη Rmax /(fm·GeV) (C5) based on the other fits.
From Tab. XI, we observe that the fits to the long dis-
• fit6: tance parts of “48I light con” and “48I light discon” are
very unreliable. Note that the statistical error for the
6
Rmax connected diagrams and disconnected diagrams are al-
f (Rmax ) = A/fm 3
e−mη Rmax /(fm·GeV)
(Rmax + (C fm)3 )2 most independent. If we fit the contributions from the
(C6) connected and disconnected diagrams individually and
 21

add the results, the statistical error will approximately

add up in quadrature. The difficulty of fitting the contri-
butions from the individual connected and disconnected
diagrams is likely due to the long tail from π 0 exchange.
The large statistical and finite volume error in the long
distance region make it hard to quantify the size of this
long tail. However, combining them to form the “light
no-pion” contribution with Eq. (29) cancels the long dis-
tance π 0 exchange contribution entirely, making the fit
much more reliable. Still, we found some fitting form de-
pendence even after the cancellation of the π 0 exchange
contribution, which we accommodate with the 100% sys-
tematic uncertainty estimated previously.
 22

Contribution name Form fit range/fm χ2d.o.f. p-value A B C

48I light no-pion fit1 0.5 - 4.0 0.94 0.305 130.57748 0.62982 0.66432
48I light no-pion fit2 1.0 - 4.0 0.85 0.481 102.07658 0.61515 -
48I light no-pion fit3 0.5 - 4.0 1.20 0.110 25.57747 0.24408 0.70365
48I light no-pion fit4 0.5 - 4.0 1.04 0.274 76.96882 - 0.55112
48I light no-pion fit5 1.0 - 4.0 0.95 0.433 69.98107 - -
48I light no-pion fit6 0.5 - 4.0 1.56 0.031 275.93065 - 0.98944
48I light no-pion fit7 0.5 - 4.0 1.18 0.160 32.61837 - 0.73081
48I strange con fit1 0.5 - 4.0 17.72 0.000 18.07733 0.76754 0.64582
48I strange con fit2 1.0 - 4.0 39.77 0.000 13.62403 0.74617 -
48I strange con fit3 0.5 - 4.0 79.03 0.000 15.87426 0.53922 0.94749
48I light con fit1 1.5 - 4.0 1.25 0.088 12.77979 0.26616 -0.59554
48I light con fit2 1.5 - 4.0 1.20 0.168 14.56237 0.27545 -
48I light con fit3 1.5 - 4.0 1.67 0.024 50.98630 0.13900 1.34340
48I light con fit8 1.5 - 4.0 1.60 0.042 50.98629 - 1.34340
48I light discon fit1 1.5 - 4.0 1.17 0.162 -43.31953 0.32417 2.85847
48I light discon fit2 1.5 - 4.0 1.48 0.060 -1.92375 0.16497 -
48I light discon fit3 1.5 - 4.0 1.12 0.212 -495.03384 0.22956 2.81122
48I light discon fit8 1.5 - 4.0 1.10 0.210 -72.47126 - 2.18553

Contribution name Form Rmax > 4fm Rmax > 2.5fm Rmax > 2.0fm
48I light no-pion fit1 0.01(0.02) 0.32(0.24) 0.89(0.48)
48I light no-pion fit2 0.01(0.02) 0.32(0.21) 0.85(0.40)
48I light no-pion fit3 0.15(0.27) 0.92(0.76) 1.67(1.02)
48I light no-pion fit4 0.03(0.00) 0.63(0.04) 1.46(0.09)
48I light no-pion fit5 0.03(0.00) 0.58(0.04) 1.34(0.08)
48I light no-pion fit6 0.00(0.00) 0.09(0.01) 0.32(0.05)
48I light no-pion fit7 0.08(0.01) 0.66(0.06) 1.31(0.12)
48I strange con ×1000 fit1 0.06(0.00) 5.90(0.06) 22.72(0.15)
48I strange con ×1000 fit2 0.08(0.00) 6.13(0.05) 22.64(0.14)
48I strange con ×1000 fit3 0.10(0.00) 5.79(0.05) 21.40(0.14)
48I light con fit1 4.96(2.11) 13.14(3.04) 16.69(3.08)
48I light con fit2 4.43(1.06) 12.40(1.86) 15.96(2.00)
48I light con fit3 4.14(0.48) 11.01(1.00) 14.54(1.17)
48I light con fit8 4.14(0.39) 11.01(0.93) 14.54(1.12)
48I light discon fit1 −3.13(4.06) −9.21(5.04) −11.20(5.07)
48I light discon fit2 −13.23(6.56) −19.61(7.61) −21.27(7.70)
48I light discon fit3 −2.76(2.33) −8.73(3.50) −10.78(3.57)
48I light discon fit8 −5.22(1.84) −11.46(3.08) −13.46(3.24)

TABLE XI. The long distance part of the ano-pion

µ × 1010 as defined in Eq. 29 from different fits to the 48I ensemble data.
The numbers in parentheses are statistical uncertainties. We perform uncorrelated fit for the summand. The minimized χ2d.o.f.
is shown in the table. We also calculated the p-value using the method invented in Ref. [74], which allows a correct p-value
determined for uncorrelated fits. In addition to the light no-pion contribution, we also perform the relevant fit to the data from
stanage quark connected diagrams, light quark connected diagrams, and light quark disconnected diagrams. The contributions
from the 48I strange quark connected diagrams are much smaller compared with other quantities. We multiply them by a
factor of 1000 to fit them in this table.
 23

48I light no-pion summand 48I light no-pion

fit1 fit1
15 12

10
10
8

aµ × 1010 partial-sum
aµ × 1010 summand 5
6

4
0

2
−5
0

−10 −2
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 12. Same as Fig. 6, reproduced for easy comparison with the other figures in this appendix. The fit function is Eq. (C1).
The fit starts at 0.5 fm.

48I light no-pion summand 48I light no-pion

fit2 fit2
15 12

10
10
8
aµ × 1010 partial-sum
aµ × 1010 summand

5
6

4
0

2
−5
0

−10 −2
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 13. Similar to Fig. 6 but with the fit function in Eq. (C2). The fit starts at 1.0 fm.

48I light no-pion summand 48I light no-pion

fit3 fit3
15 12

10
10
8
aµ × 1010 partial-sum
aµ × 1010 summand

5
6

4
0

2
−5
0

−10 −2
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 14. Similar to Fig. 6 but with the fit function in Eq. (C3). The fit starts at 0.5 fm.
 24

48I light no-pion summand 48I light no-pion

fit4 fit4
15 12

10
10
8

aµ × 1010 partial-sum
aµ × 1010 summand 5
6

4
0

2
−5
0

−10 −2
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 15. Similar to Fig. 6 but with the fit function in Eq. (C4). The fit starts at 0.5 fm.

48I light no-pion summand 48I light no-pion

fit5 fit5
15 12

10
10
8
aµ × 1010 partial-sum
aµ × 1010 summand

5
6

4
0

2
−5
0

−10 −2
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 16. Similar to Fig. 6 but with the fit function in Eq. (C5). The fit starts at 1.0 fm.

48I light no-pion summand 48I light no-pion

fit6 fit6
15 12

10
10
8
aµ × 1010 partial-sum
aµ × 1010 summand

5
6

4
0

2
−5
0

−10 −2
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 17. Similar to Fig. 6 but with the fit function in Eq. (C6). The fit starts at 0.5 fm.
 25

48I light no-pion summand 48I light no-pion

fit7 fit7
15 12

10
10
8

aµ × 1010 partial-sum
aµ × 1010 summand
5
6

4
0

2
−5
0

−10 −2
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 18. Similar to Fig. 6 but with the fit function in Eq. (C7). The fit starts at 0.5 fm.

48I strange con summand 48I strange con

fit1 fit1
0.3 0.35

0.25 0.3

0.25
0.2
aµ × 1010 partial-sum
aµ × 1010 summand

0.2
0.15
0.15
0.1
0.1
0.05
0.05

0 0

−0.05 −0.05
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 19. Similar to Fig. 6 but for the strange quark connected contriubtion with the fit function in Eq. (C1). The fit starts at
0.5 fm.
48I strange con summand 48I strange con
fit2 fit2
0.35 0.35

0.3 0.3

0.25 0.25
aµ × 1010 partial-sum
aµ × 1010 summand

0.2 0.2

0.15 0.15

0.1 0.1

0.05 0.05

0 0

−0.05 −0.05
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 20. Similar to Fig. 6 but for the strange quark connected contribution with the fit function in Eq. (C2). The fit starts at
1.0 fm.
 26

48I strange con summand 48I strange con

fit3 fit3
0.35 0.35

0.3 0.3

0.25 0.25

aµ × 1010 partial-sum
aµ × 1010 summand
0.2 0.2

0.15 0.15

0.1 0.1

0.05 0.05

0 0

−0.05 −0.05
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 21. Similar to Fig. 6 but for the strange quark connected contribution with the fit function in Eq. (C3). The fit starts at
0.5 fm.
48I light con summand 48I light con
fit1 fit1
12 25

10
20
8
aµ × 1010 partial-sum
aµ × 1010 summand

6 15

4
10
2

0 5

−2
0
−4

−6 −5
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 22. Similar to Fig. 6 but for light quark connected contribution with the fit function in Eq. (C1). The fit starts at 1.5 fm.

48I light con summand 48I light con

fit2 fit2
12 25

10
20
8
aµ × 1010 partial-sum
aµ × 1010 summand

6 15

4
10
2

0 5

−2
0
−4

−6 −5
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 23. Similar to Fig. 6 but for light quark connected contribution with a the fit function in Eq. (C2). The fit starts at
1.5 fm.
 27

48I light con summand 48I light con

fit3 fit3
12 25

10
20
8

aµ × 1010 partial-sum
aµ × 1010 summand

6 15

4
10
2

0 5

−2
0
−4

−6 −5
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 24. Similar to Fig. 6 but for the light quark connected contribution with the fit function in Eq. (C3). The fit starts at
1.5 fm.

48I light con summand 48I light con

fit8 fit8
12 25

10
20
8
aµ × 1010 partial-sum
aµ × 1010 summand

6 15

4
10
2

0 5

−2
0
−4

−6 −5
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 25. Similar to Fig. 6 but for the light quark connected contribution with the fit function in Eq. (C8). The fit starts at
1.5 fm.
 28

48I light discon summand 48I light discon

fit1 fit1
15 2

0
10
−2

aµ × 1010 partial-sum
aµ × 1010 summand

5
−4

0 −6

−8
−5
−10
−10
−12

−15 −14
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 26. Similar to Fig. 6 but for the light quark disconnected contribution with the fit function in Eq. (C1). The fit starts at
1.5 fm.

48I light discon summand 48I light discon

fit2 fit2
15 2

0
10
−2
aµ × 1010 partial-sum
aµ × 1010 summand

5 −4

−6
0
−8

−5 −10

−12
−10
−14

−15 −16
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 27. Similar to Fig. 6 but for the light quark disconnected conrtibution with the fit function in Eq. (C2). The fit starts at
1.5 fm.
 29

48I light discon summand 48I light discon

fit3 fit3
15 2

0
10
−2

aµ × 1010 partial-sum
aµ × 1010 summand

5
−4

0 −6

−8
−5
−10
−10
−12

−15 −14
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 28. Similar to Fig. 6 but for the light quark disconnected contribution with the fit function in Eq. (C3). The fit starts at
1.5 fm.

48I light discon summand 48I light discon

fit8 fit8
15 2

0
10
−2
aµ × 1010 partial-sum
aµ × 1010 summand

5
−4

0 −6

−8
−5
−10
−10
−12

−15 −14
0 1 2 3 4 5 0 1 2 3 4 5
Rmax (fm) Rmax (fm)

FIG. 29. Similar to Fig. 6 but for the light quark disconnected contribution with the fit function in Eq. (C8). The fit starts at
1.5 fm.
 30

[1] B. Abi et al. (Muon g-2), L. Laub, and P. Stoffer, JHEP 03, 101 (2020),
Phys. Rev. Lett. 126, 141801 (2021), arXiv:1910.13432 [hep-ph].
arXiv:2104.03281 [hep-ex]. [24] T. Blum, N. Christ, M. Hayakawa,
[2] M. Abe et al., PTEP 2019, 053C02 (2019), T. Izubuchi, L. Jin, C. Jung, and
arXiv:1901.03047 [physics.ins-det]. C. Lehner, Phys. Rev. Lett. 124, 132002 (2020),
[3] Y. Sato (J-PARC E34), arXiv:1911.08123 [hep-lat].
JPS Conf. Proc. 33, 011110 (2021). [25] G. Colangelo, M. Hoferichter, A. Nyffeler, M. Passera,
[4] G. Bennett et al. (Muon G- and P. Stoffer, Phys. Lett. B735, 90 (2014),
2), Phys.Rev. D73, 072003 (2006), arXiv:1403.7512 [hep-ph].
arXiv:hep-ex/0602035 [hep-ex]. [26] J. Prades, E. de Rafael, and A. Vainshtein,
[5] T. Aoyama et al., Phys. Rept. 887, 1 (2020), Adv. Ser. Direct. High Energy Phys. 20, 303 (2009),
arXiv:2006.04822 [hep-ph]. arXiv:0901.0306 [hep-ph].
[6] T. Aoyama, M. Hayakawa, T. Kinoshita, and [27] A. Nyffeler, Phys. Rev. D 79, 073012 (2009),
M. Nio, Phys. Rev. Lett. 109, 111808 (2012), arXiv:0901.1172 [hep-ph].
arXiv:1205.5370 [hep-ph]. [28] F. Jegerlehner and A. Nyffeler,
[7] T. Aoyama, T. Kinoshita, and M. Nio, Phys. Rept. 477, 1 (2009), arXiv:0902.3360 [hep-ph].
Atoms 7, 28 (2019). [29] F. Jegerlehner, EPJ Web Conf. 166, 00022 (2018),
[8] A. Czarnecki, W. J. Marciano, and A. Vainshtein, arXiv:1705.00263 [hep-ph].
Phys. Rev. D67, 073006 (2003), [Erratum: Phys. Rev. [30] J. Green, O. Gryniuk, G. von Hippel, H. B. Meyer,
D73, 119901 (2006)], arXiv:hep-ph/0212229 [hep-ph]. and V. Pascalutsa, Phys. Rev. Lett. 115, 222003 (2015),
[9] C. Gnendiger, D. Stöckinger, and arXiv:1507.01577 [hep-lat].
H. Stöckinger-Kim, Phys. Rev. D88, 053005 (2013), [31] A. Gérardin, J. Green, O. Gryniuk, G. von
arXiv:1306.5546 [hep-ph]. Hippel, H. B. Meyer, V. Pascalutsa, and
[10] M. Davier, A. Hoecker, B. Malaescu, and H. Wittig, Phys. Rev. D 98, 074501 (2018),
Z. Zhang, Eur. Phys. J. C77, 827 (2017), arXiv:1712.00421 [hep-lat].
arXiv:1706.09436 [hep-ph]. [32] V. Pauk and M. Vanderhaeghen,
[11] A. Keshavarzi, D. Nomura, and Eur. Phys. J. C 74, 3008 (2014),
T. Teubner, Phys. Rev. D97, 114025 (2018), arXiv:1401.0832 [hep-ph].
arXiv:1802.02995 [hep-ph]. [33] I. Danilkin and M. Vander-
[12] G. Colangelo, M. Hoferichter, and P. Stoffer, haeghen, Phys. Rev. D 95, 014019 (2017),
JHEP 02, 006 (2019), arXiv:1810.00007 [hep-ph]. arXiv:1611.04646 [hep-ph].
[13] M. Hoferichter, B.-L. Hoid, and B. Kubis, [34] F. Jegerlehner, The Anomalous Magnetic Moment of the Muon,
JHEP 08, 137 (2019), arXiv:1907.01556 [hep-ph]. Vol. 274 (Springer, Cham, 2017).
[14] M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, [35] M. Knecht, S. Narison, A. Rabemananjara, and
Eur. Phys. J. C80, 241 (2020), [Erratum: Eur. Phys. J. D. Rabetiarivony, Phys. Lett. B 787, 111 (2018),
C80, 410 (2020)], arXiv:1908.00921 [hep-ph]. arXiv:1808.03848 [hep-ph].
[15] A. Keshavarzi, D. Nomura, and T. Teub- [36] G. Eichmann, C. S. Fischer, and
ner, Phys. Rev. D101, 014029 (2020), R. Williams, Phys. Rev. D 101, 054015 (2020),
arXiv:1911.00367 [hep-ph]. arXiv:1910.06795 [hep-ph].
[16] A. Kurz, T. Liu, P. Marquard, and [37] P. Roig and P. Sanchez-Puertas,
M. Steinhauser, Phys. Lett. B734, 144 (2014), Phys. Rev. D 101, 074019 (2020),
arXiv:1403.6400 [hep-ph]. arXiv:1910.02881 [hep-ph].
[17] K. Melnikov and A. Vain- [38] A. Gérardin, H. B. Meyer, and A. Nyf-
shtein, Phys. Rev. D70, 113006 (2004), feler, Phys. Rev. D 94, 074507 (2016),
arXiv:hep-ph/0312226 [hep-ph]. arXiv:1607.08174 [hep-lat].
[18] P. Masjuan and P. Sanchez- [39] C. Alexandrou et al., (2022), arXiv:2212.06704 [hep-lat].
Puertas, Phys. Rev. D 95, 054026 (2017), [40] S. Burri et al., (2022), arXiv:2212.10300 [hep-lat].
arXiv:1701.05829 [hep-ph]. [41] T. Blum, S. Chowdhury, M. Hayakawa, and T. Izubuchi,
[19] G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer, (2014), arXiv:1407.2923 [hep-lat].
JHEP 04, 161 (2017), arXiv:1702.07347 [hep-ph]. [42] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi,
[20] M. Hoferichter, B.-L. Hoid, B. Kubis, S. Le- L. Jin, and C. Lehner, Phys. Rev. D 93, 014503 (2016),
upold, and S. P. Schneider, JHEP 10, 141 (2018), arXiv:1510.07100 [hep-lat].
arXiv:1808.04823 [hep-ph]. [43] T. Blum, N. Christ, M. Hayakawa,
[21] A. Gérardin, H. B. Meyer, and A. Nyf- T. Izubuchi, L. Jin, C. Jung, and
feler, Phys. Rev. D 100, 034520 (2019), C. Lehner, Phys. Rev. Lett. 118, 022005 (2017),
arXiv:1903.09471 [hep-lat]. arXiv:1610.04603 [hep-lat].
[22] J. Bijnens, N. Hermansson- [44] M. Hayakawa and S. Uno,
Truedsson, and A. Rodríguez- Prog. Theor. Phys. 120, 413 (2008),
Sánchez, Phys. Lett. B798, 134994 (2019), arXiv:0804.2044 [hep-ph].
arXiv:1908.03331 [hep-ph]. [45] N. Asmussen, J. Green, H. B. Meyer, and
[23] G. Colangelo, F. Hagelstein, M. Hoferichter, A. Nyffeler, PoS LATTICE2016, 164 (2016),
 31

arXiv:1609.08454 [hep-lat]. arXiv:1908.07050 [hep-lat].

[46] T. Blum, N. Christ, M. Hayakawa, [60] Y. Li, S.-C. Xia, X. Feng, L.-C. Jin, and
T. Izubuchi, L. Jin, C. Jung, and C. Liu, Phys. Rev. D 103, 014514 (2021),
C. Lehner, Phys. Rev. D 96, 034515 (2017), arXiv:2009.01029 [hep-lat].
arXiv:1705.01067 [hep-lat]. [61] J. Bijnens and J. Relefors, JHEP 09, 113 (2016),
[47] E.-H. Chao, A. Gérardin, J. R. Green, R. J. Hudsp- arXiv:1608.01454 [hep-ph].
ith, and H. B. Meyer, Eur. Phys. J. C 80, 869 (2020), [62] L. Jin, T. Blum, N. Christ, M. Hayakawa,
arXiv:2006.16224 [hep-lat]. T. Izubuchi, C. Jung, and C. Lehner,
[48] E.-H. Chao, R. J. Hudspith, A. Gérardin, PoS LATTICE2016, 181 (2016),
J. R. Green, H. B. Meyer, and arXiv:1611.08685 [hep-lat].
K. Ottnad, Eur. Phys. J. C 81, 651 (2021), [63] One can use different setup for the connected and dis-
arXiv:2104.02632 [hep-lat]. connected diagrams, in which case this relation will not
[49] E.-H. Chao, R. J. Hudspith, A. Gérardin, J. R. Green, hold exactly. For example, Ref. [48] uses different setups.
and H. B. Meyer, Eur. Phys. J. C 82, 664 (2022), [64] Y. Zhao, Lattice Calculation of the π 0 → e+ e− and the
arXiv:2204.08844 [hep-lat]. KL → γγ Decays, Ph.D. thesis, Columbia U. (2022).
[50] S. Laporta and E. Remiddi, [65] N. Asmussen, E.-H. Chao, A. Gérardin, J. R. Green,
Phys. Lett. B 301, 440 (1993). R. J. Hudspith, H. B. Meyer, and A. Nyffeler, (2022),
[51] S. Laporta and E. Remiddi, arXiv:2210.12263 [hep-lat].
Phys. Lett. B 265, 182 (1991). [66] X. Feng, S. Aoki, H. Fukaya, S. Hashimoto,
[52] T. Blum et al. (RBC, UKQCD), T. Kaneko, J.-i. Noaki, and E. Shin-
Phys. Rev. D 93, 074505 (2016), tani, Phys. Rev. Lett. 109, 182001 (2012),
arXiv:1411.7017 [hep-lat]. arXiv:1206.1375 [hep-lat].
[53] D. Renfrew, T. Blum, N. Christ, R. Mawhin- [67] B. Moussallam, Phys. Rev. D 51, 4939 (1995),
ney, and P. Vranas, PoS LATTICE2008, 048 (2008), arXiv:hep-ph/9407402.
arXiv:0902.2587 [hep-lat]. [68] M. Knecht, S. Peris, M. Perrottet, and
[54] M. Abramczyk, S. Aoki, T. Blum, T. Izubuchi, H. Ohki, E. de Rafael, Phys. Rev. Lett. 83, 5230 (1999),
and S. Syritsyn, Phys. Rev. D 96, 014501 (2017). arXiv:hep-ph/9908283.
[55] T. Blum, P. A. Boyle, V. Gülpers, T. Izubuchi, [69] M. J. Ramsey-Musolf and M. B.
L. Jin, C. Jung, A. Jüttner, C. Lehner, Wise, Phys. Rev. Lett. 89, 041601 (2002),
A. Portelli, and J. T. Tsang (RBC, arXiv:hep-ph/0201297.
UKQCD), Phys. Rev. Lett. 121, 022003 (2018), [70] T. Kinoshita, B. Nizic, and Y. Okamoto,
arXiv:1801.07224 [hep-lat]. Phys. Rev. D 31, 2108 (1985).
[56] T. Blum et al., (2023), arXiv:2301.08696 [hep-lat]. [71] M. Hayakawa, T. Kinoshita, and A. I. Sanda,
[57] M. A. Clark, C. Jung, and C. Lehner, in Phys. Rev. Lett. 75, 790 (1995), arXiv:hep-ph/9503463.
[72]2017)
35th International Symposium on Lattice Field Theory (Lattice J. Granada,
Bijnens,Spain,EPJ
JuneWeb Conf.
18-24, 118, 01002 (2016),
2017
(2017) arXiv:1710.06884 [hep-lat]. arXiv:1510.05796 [hep-ph].
[58] T. Blum, T. Izubuchi, and E. Shin- [73] N. Christ, X. Feng, L. Jin, C. Tu, and Y. Zhao, (2022),
tani, Phys. Rev. D 88, 094503 (2013), arXiv:2208.03834 [hep-lat].
arXiv:1208.4349 [hep-lat]. [74] N. Christ, R. Eranki, and C. Kelly (RBC, UKQCD),
[59] W. Detmold, D. J. Murphy, A. V. Pochin- (2024), arXiv:2409.11379 [hep-lat].
sky, M. J. Savage, P. E. Shanahan, and
M. L. Wagman, Phys. Rev. D 104, 034502 (2021),