RBC and UKQCD Collaborations

#  Hadronic light-by-light contribution to the muon anomaly from lattice QCD
with infinite volume QED at physical pion mass

Thomas Blum<EMAIL_ADDRESS>Physics Department, University of
Connecticut, Storrs, CT 06269-3046, USA  Norman Christ  Physics Department,
Columbia University, New York, NY 10027, USA  Masashi Hayakawa  Department of
Physics, Nagoya University, Nagoya 464-8602, Japan  Nishina Center, RIKEN,
Wako, Saitama 351-0198, Japan  Taku Izubuchi  Physics Department, Brookhaven
National Laboratory, Upton, NY 11973, USA  RIKEN-BNL Research Center,
Brookhaven National Laboratory, Upton, NY 11973, USA  Luchang Jin
<EMAIL_ADDRESS>Physics Department, University of Connecticut, Storrs,
CT 06269-3046, USA  RIKEN-BNL Research Center, Brookhaven National Laboratory,
Upton, NY 11973, USA  Chulwoo Jung  Physics Department, Brookhaven National
Laboratory, Upton, NY 11973, USA  Christoph Lehner  Fakultät für Physik,
Universität Regensburg, Universitätsstraße 31, 93040 Regensburg, Germany
Cheng Tu  Physics Department, University of Connecticut, Storrs, CT
06269-3046, USA

###### Abstract

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 $a_{\mu}^{\text{HLbL}}=12.47(1.15)(0.99)\times
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~{}\mathrm{GeV}$, lattice size
$L=5.5~{}\mathrm{fm}$, and $m_{\pi}=139~{}\mathrm{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

Muons are spin-$1/2$ charged particles with non-zero magnetic moment:

$\displaystyle\boldsymbol{\mu}=-g\frac{e}{2m}\bf S,$ (1)

where $\bf S$ is the particle’s spin, $e$ and $m$ are the electric charge and
mass, respectively, and $g$ is the Landé $g$-factor. The Dirac equation
predicts that $g=2$, exactly, so any difference from 2 must arise from
interactions. The magnetic moment of a fermion can be defined in terms of its
electromagnetic form factors in the zero momentum transfer limit. Lorentz and
gauge symmetries tightly constrain the form of the interactions. In Euclidean
space-time:

$\displaystyle\langle\mu({\bf p^{\prime}},\,s^{\prime})|J_{\nu}(0)|\mu({\bf
p},\,s)\rangle$ (2) $\displaystyle=-e\bar{u}_{s^{\prime}}({\bf
p^{\prime}})\left(F_{1}(q^{2})\gamma_{\nu}+i\frac{F_{2}(q^{2})}{4m}[\gamma_{\nu},\gamma_{\rho}]q_{\rho}\right)u_{s}({\bf
p}),$

where $J_{\nu}$ is the electromagnetic current, and $F_{1}$ and $F_{2}$ are
form factors, giving the charge and magnetic moment at zero momentum transfer
($q^{2}=(p^{\prime}-p)^{2}=0$), or static limit. The ${u}_{s}(\bf{p})$ and
$\bar{u}_{s}({\bf p})$ are Dirac spinors. The anomalous part of the magnetic
moment is given by $F_{2}(0)$ alone, and is known as the anomaly,

$\displaystyle a_{\mu}\equiv F_{2}(q^{2}=0)=\frac{g-2}{2}.$ (3)

The muon anomalous magnetic moment is one of the most precisely determined
quantities in particle physics. Compared with the electron anomalous magnetic
moment, which is determined with higher accuracy, the muon is expected to be
much more sensitive to new physics at very large energy scales due to its
heavier mass. Two experiments Fermilab (E989) [1] and J-PARC (E34) [2, 3] are
aiming at even higher precision in measuring the muon $g-2$. The initial
result released by Fermilab (E989) [1] confirmed the previously best result
obtained by the BNL E821 experiment [4] and reduced the experimental
uncertainty from 0.54 ppm to 0.46 ppm. The final goal of the Fermilab
experiment is to reduce the uncertainty further to approximately 0.14 ppm. The
J-PARC experiment adopts a very different measurement strategy and will serve
as an important cross-check. Its final accuracy goal is 0.45 ppm for
statistical uncertainty and 0.07 ppm for systematic uncertainty.

The Standard Model result provided by the Muon $g-2$ Theory Initiative [5, 6,
7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]
currently has an uncertainty of 0.37 ppm and is in $4.2\sigma$ tension with
the experimental value. The two leading sources of uncertainty come from the
leading QCD contributions, hadronic vacuum polarization (HVP) and hadronic
light-by-light (HLbL) scattering, both illustrated in Fig. 1.

Figure 1:  The HVP (left) and HLbL (right) contribution to muon $g-2$. The
shaded circles represent hadronic interactions.

For a long time, the HLbL contribution has been estimated by hadronic models
[26, 27, 28, 29]. In Ref. [30, 31], the hadronic light-by-light scattering
amplitudes estimated by hadronic models are compared with lattice QCD
calculations. The result from hadronic models is difficult to improve further.
More recently, the HLbL contribution has been calculated with the dispersion
relation approach [17, 18, 19, 20, 21, 22, 23, 32, 33, 34, 35, 36, 37, 24,
25]. These results are summarized in the 2020 white paper [5],

$\displaystyle a_{\mu}^{\text{HLbL,WP2020,pheno}}\times 10^{10}=9.2(1.9)$ (4)

Within the dispersion relation framework, the pseudo-scalar pole is the
leading contribution to HLbL. This contribution is defined in terms of the
pseudo-scalar transition form factors to two (possibly off-shell) photons,
$\pi^{0},\eta,\eta^{\prime}\to\gamma^{*}\gamma^{*}$. These form factors can
also be calculated with lattice QCD. [38, 21, 39, 40]

The first direct lattice QCD calculation of the HLbL contribution was
performed by the RBC-UKQCD collaborations. [41] Then, we made many
improvements to the calculation method and also included the contribution from
the leading disconnected diagrams. [42, 43] In our previous work [24], we
applied the method developed to lattice gauge ensembles with different lattice
spacings and spatial volumes, all at physical pion mass. The calculations
incorporate both QCD and QED in the same finite volume lattice using the QEDL
scheme [44]. We extrapolated the results to infinite volume and the continuum
limit. The final result was

$\displaystyle a_{\mu}^{\text{HLbL,RBC-UKQCD2019}}\times
10^{10}=7.87(3.06)_{\text{stat}}(1.77)_{\text{syst}}\,.$ (5)

This was the first lattice QCD result for HLbL with all systematic effects
controlled. The statistical error was still larger than the phenomenological
results obtained with dispersion relations or the previous hadronic model
estimation.

The Mainz group pioneered the QED∞ approach and pre-calculated the QED kernel
semi-analytically in infinite volume [30, 45]. We also explored this strategy
and independently calculated the infinite volume QED weighting function [46].
In that work, we found that while the QED∞ approach ensures an exponentially
suppressed finite volume error, the size of the finite volume error can still
be significant. However, making use of the current conservation property of
the hadronic four point function, we designed subtractions to the infinite
volume QED weighting functions to reduce the finite volume errors and also the
discretization errors. The subtracted infinite volume QED weighting function
is used in the current work.

The effectiveness of the subtraction is confirmed by the Mainz group and a
different subtraction scheme was adopted in their calculation. [47, 48, 49,
45]. Their calculation was performed with pion masses ranging from
$200~{}\mathrm{MeV}$ to $420~{}\mathrm{MeV}$ and extrapolated to the physical
point. All sub-leading disconnected diagrams were carefully calculated and
were found to be consistent with zero. The charm quark contribution from the
connected and disconnected diagrams were calculated as well. Finally, a
statistically more precise result was obtained:

$\displaystyle a_{\mu}^{\text{HLbL,Mainz2021}}\times 10^{10}=10.96(1.59).$ (6)

In this work, we present our latest lattice calculation of the HLbL
contribution to muon $g-2$ with subtracted infinite volume QED weighting
function as developed in our previous work [46]. The calculation is directly
performed at physical pion mass and therefore eliminates the systematic
uncertainty from chiral extrapolations. We find that calculating the HLbL
contribution directly at the physical pion mass is considerably more difficult
than a calculation at a heavier pion mass due the larger contribution and
statistical fluctuations from the long distance region. Following the
suggestion of Ref. [48], we focus on the calculation of the connected and
leading disconnected diagrams. The main calculation is performed with a single
lattice spacing where $a^{-1}=1.73~{}\mathrm{GeV}$. While we use Möbius domain
wall fermions (MDWF) to eliminate $\mathcal{O}(a)$ discretization errors, the
remaining $\mathcal{O}(a^{2})$ effects are estimated to be the largest source
of systematic uncertainty of this calculation.

The paper is organized as follows: in section II and III, we review the
calculation method developed in our previous works and describe new techniques
used in this calculation. We then show our numerical results in section IV and
conclude in section V.

##  II Theoretical framework and calculation method

We base our current calculation on the framework already set up in our
previous works [42, 43, 46]. In this work, we combine the method used in Ref.
[42, 43] to calculate the hadronic four point function with lattice QCD with
the subtracted infinite volume QED weighting function developed in Ref. [46].
We have tested the QED weighting function by calculating the leptonic light-
by-light contribution to muon $g-2$ and have reproduced the results obtained
by analytical calculation [50, 51]. We outline the method below and focus on
the improvements we made in this work.

Figure 2: Diagrams contributing to the muon anomaly. The diagrams in the top
row are the leading ones, and do not vanish in the $SU(3)$ flavor limit.
Strong interactions to all orders, including gluons connecting the quark loops
and sea quark loops which are not connected by photons, are not shown.

All of the diagrams contributing at $O(\alpha^{3})$ to HLbL scattering are
shown in Fig. 2. We use the term (quark-) “connected diagram” to refer to the
one on the left on the top row and “disconnected diagram” to refer to the one
on the right. The remaining diagrams, which are all suppressed in the flavor
SU(3) limit, are referred to as the “sub-leading disconnected diagrams”. We
only explicitly draw quark loops that are connected to photons. Gluons and sea
quark loops that are not connected to photons are not shown in the figure but
are included automatically in dynamical lattice QCD calculations. To make non-
zero contributions to the muon $g-2$, the quark loops in the disconnected
diagrams must be connected by gluons.

The contribution to the muon $g-2$ can be calculated with the combination of
the hadronic four-point function $\mathcal{H}$ and the QED weighting function
$\mathcal{G}$ [46]:

$\displaystyle
a_{\mu}^{\text{HLbL}}\frac{e}{m}\bar{u}_{s^{\prime}}(\vec{0})\frac{\Sigma_{i}}{2}u_{s}(\vec{0})$
(7)
$\displaystyle=\frac{1}{VT}\sum_{x_{\text{op}}}\sum_{x,y,z}\frac{1}{2}\epsilon_{i,j,k}\big{(}x_{\text{op}}-x_{\text{ref}}(x,y,z)\big{)}_{j}$
$\displaystyle\hskip 5.69046pt\times
i^{3}e^{6}\mathcal{H}_{k,\rho,\sigma,\lambda}(x_{\text{op}},x,y,z)\bar{u}_{s^{\prime}}(\vec{0})\mathcal{G}_{\rho,\sigma,\kappa}(x,y,z)u_{s}(\vec{0}),$

where $\bar{u}_{s^{\prime}}(\vec{0})$, $u_{s}(\vec{0})$ are Dirac spinors for
the outgoing and incoming muon in the diagram, respectively.
$\Sigma_{k}=\epsilon_{i,j,k}\gamma_{i}\gamma_{j}/(2i)$ is the $4\times 4$
version of the Pauli matrix, $\sigma_{k}$. From the spin structure of the muon
particle, we can obtain the expression for $a_{\mu}^{\text{HLbL}}$:

$\displaystyle a_{\mu}^{\text{HLbL}}$
$\displaystyle=\frac{2me^{2}}{3}\frac{1}{VT}\sum_{x_{\text{op}}}\sum_{x,y,z}\frac{1}{2}\epsilon_{i,j,k}\big{(}x_{\text{op}}-x_{\text{ref}}(x,y,z)\big{)}_{j}$
$\displaystyle\times
6e^{4}\mathcal{H}_{k,\rho,\sigma,\lambda}(x_{\text{op}},x,y,z)\mathcal{M}_{i,\rho,\sigma,\lambda}(x,y,z)$
(8)

where

$\displaystyle\mathcal{M}_{i,\rho,\sigma,\lambda}(x,y,z)=\frac{1}{2}\mathrm{Tr}\Big{[}\frac{1}{6}i^{3}\mathcal{G}_{\rho,\sigma,\nu}(x,y,z)\Sigma_{i}\Big{]}\,.$
(9)

Figure 3: Diagramatic representation of the QED weighting function defined in
Eq. 9, following Ref. [46].

The QED weighting function $\mathcal{G}$ is shown diagramatically in Fig. 3
and is expressed in terms of the free muon and Feynman gauge photon
propagators, $S_{{\mu}}$ and $G$:

$\displaystyle\mathfrak{G}_{\sigma,\kappa,\rho}(y,z,x)$
$\displaystyle=\lim_{t_{\text{src}}\rightarrow-\infty,t_{\text{snk}}\rightarrow\infty}e^{m_{{\mu}}\left(t_{\text{snk}}-t_{\text{src}}\right)}$
(10)
$\displaystyle\times\int_{\alpha,\beta,\eta,\vec{x}_{\text{snk}},\vec{x}_{\text{src}}}G(x,\alpha)G(y,\beta)G(z,\eta)$
$\displaystyle\times
S_{{\mu}}\left(x_{\text{snk}},\beta\right)i\gamma_{\sigma}S_{{\mu}}(\beta,\eta)i\gamma_{\kappa}S_{{\mu}}(\eta,\alpha)i\gamma_{\rho}S_{{\mu}}\left(\alpha,x_{\text{src}}\right)\,.$

As is well known, the above expression contains an infrared divergence that
vanishes after projection to its magnetic part. We can also remove this
infrared divergent piece by the following procedure:

$\displaystyle\mathfrak{G}^{(1)}_{\sigma,\kappa,\rho}(y,z,x)$
$\displaystyle=\frac{1}{2}\mathfrak{G}_{\sigma,\kappa,\rho}(y,z,x)$ (11)
$\displaystyle+\frac{1}{2}[\mathfrak{G}_{\rho,\kappa,\sigma}(x,z,y)]^{\dagger}\,.$

In addition we can perform somewhat arbitrary subtractions to this infinite
volume QED weighting function without changing the final result due to vector
current conservation satisfied by the hadronic four point function,

$\displaystyle\mathfrak{G}^{(2)}_{\sigma,\kappa,\rho}(y,z,x)$
$\displaystyle=\mathfrak{G}^{(1)}_{\sigma,\kappa,\rho}(y,z,x)$ (12)
$\displaystyle-\mathfrak{G}^{(1)}_{\sigma,\kappa,\rho}(z,z,x)-\mathfrak{G}^{(1)}_{\sigma,\kappa,\rho}(y,z,z).$

Note that $\mathfrak{G}^{(1)}_{\sigma,\kappa,\rho}(z,z,z)=0$, so this
subtraction significantly reduces the size of the QED weighting function when
$|x-z|$ or $|y-z|$ is small. This is the region where the hadronic function
from the lattice calculation has the largest discretization error. It turns
out that this subtraction greatly reduces the discretization error. This is
the major finding of Ref. [46]. We should also note that the subtraction does
impact the integrand and partial sum. In Ref. [48], a modified subtraction
scheme is used, so their integrand and partial sum cannot be directly compared
with ours. Finally, we include all possible permutations of the subtracted QED
weighting function which are required for the total contribution to the muon
$g-2$:

$\displaystyle i^{3}\mathcal{G}_{\rho,\sigma,\kappa}(x,y,z)$
$\displaystyle=\mathfrak{G}^{(2)}_{\rho,\sigma,\kappa}(x,y,z)+\mathfrak{G}^{(2)}_{\sigma,\kappa,\rho}(y,z,x)$
(13)
$\displaystyle+\mathfrak{G}^{(2)}_{\kappa,\rho,\sigma}(z,x,y)+\mathfrak{G}^{(2)}_{\kappa,\sigma,\rho}(z,y,x)$
$\displaystyle+\mathfrak{G}^{(2)}_{\rho,\kappa,\sigma}(x,z,y)+\mathfrak{G}^{(2)}_{\sigma,\rho,\kappa}(y,x,z).$

Another component of the master formula Eq. II is $x_{\text{ref}}$, the
reference position for the moment method to calculate the magnetic moment.
Again, due to current conservation, there are many possible choices for
$x_{\text{ref}}$. In this work, we use the following choice for the connected
diagram:

$\displaystyle x_{\text{ref}}(x,y,z)$ $\displaystyle=x_{\text{ref-
far}}(x,y,z)$ (14) $\displaystyle=\left\\{\begin{array}[]{l}x\qquad\text{if
}|y-z|<\min(|x-y|,|x-z|)\\\ y\qquad\text{if }|x-z|<\min(|x-y|,|y-z|)\\\
z\qquad\text{if }|x-y|<\min(|x-z|,|y-z|)\\\
\frac{1}{3}(x+y+z)\qquad\text{otherwise}\end{array}\right.$ (19)

We make a slightly different choice of $x_{\text{ref}}$, Eq. (26), for the
disconnected diagrams. The rationale will be described in the later part of
this section. We use $\mathcal{H}$ to denote the hadronic four point function:

$\displaystyle 6e^{4}\mathcal{H}_{k,\rho,\sigma,\lambda}(x_{\text{op}},x,y,z)$
$\displaystyle\hskip 56.9055pt=\langle
TJ_{k}(x_{\text{op}})J_{\rho}(x)J_{\sigma}(y)J_{\lambda}(z)\rangle_{\text{QCD}}$
(20) $\displaystyle
J_{\nu}(x)=\sum_{q=u,d,s,c}e_{q}Z_{V}\bar{\psi}_{q}(x)\gamma_{\nu}\psi_{q}(x)$
(21)

where $Z_{V}$ is the lattice local vector current renormalization constant.
After Wick contraction, $\mathcal{H}$ can be expressed as the sum of different
types of diagrams as illustrated in Fig. 2. For the disconnected and sub-
leading disconnected diagrams, we require the quark loops to be connected by
gluons.

$\displaystyle\mathcal{H}=\mathcal{H}^{\text{con}}+\mathcal{H}^{\text{discon}}+\mathcal{H}^{\text{sub-
leading-discon}}+\mathcal{H}^{\text{charm}}$ (22)

where $\mathcal{H}^{\text{con}}$ includes the light and strange quark
connected diagrams, $\mathcal{H}^{\text{discon}}$ the light and strange quark
disconnected diagrams, $\mathcal{H}^{\text{sub-leading-discon}}$ the light and
strange quark sub-leading disconnected diagrams (vanish in the flavor SU(3)
limit), and $\mathcal{H}^{\text{charm}}$ all diagrams involving charm quark
loops.

Naturally, after the Wick contraction, $\mathcal{H}$ defined in Eq. (II) can
be expressed in terms of quark propagators and includes all permutations of
$x$, $y$, and $z$. However, note that all other terms in the master formula
Eq. (II) are symmetric with respect to permutations of $x$, $y$, and $z$
(along with its Lorentz indices). Therefore, we are allowed to calculate only
a subset of the diagrams generated by the Wick contractions in Eq. (II) and
multiply them with appropriate factors.

$\displaystyle\hskip
14.22636pt6e^{4}\mathcal{H}^{\text{con}}_{k,\rho,\sigma,\lambda}(x_{\text{op}},x,y,z)$
$\displaystyle\Rightarrow 6e^{4}\mathcal{H}^{\text{con-no-
perm}}_{k,\rho,\sigma,\lambda}(x_{\text{op}},x,y,z)$ (23)
$\displaystyle=-6\Big{\langle}\mathrm{Re}\sum_{q=u,d,s}e_{q}^{4}\mathrm{Tr}\Big{(}\gamma_{\nu}S_{q}(x_{\text{op}},x)\gamma_{\rho}S_{q}(x,z)$
$\displaystyle\hskip
56.9055pt\times\gamma_{\kappa}S_{q}(z,y)\gamma_{\sigma}S_{q}(y,x_{\text{op}})\Big{)}\Big{\rangle}_{\text{QCD}},$
$\displaystyle\hskip
14.22636pt6e^{4}\mathcal{H}^{\text{discon}}_{k,\rho,\sigma,\lambda}(x_{\text{op}},x,y,z)$
$\displaystyle\Rightarrow 6e^{4}\mathcal{H}^{\text{discon-no-
perm}}_{k,\rho,\sigma,\lambda}(x_{\text{op}},x,y,z)$ (24)
$\displaystyle=3\Big{\langle}\sum_{q^{\prime}=u,d,s}e_{q^{\prime}}^{2}\mathrm{Tr}\Big{(}\gamma_{\nu}S_{q^{\prime}}(x_{\text{op}},x)\gamma_{\rho}S_{q^{\prime}}(x,x_{\text{op}})\Big{)}$
$\displaystyle\hskip
28.45274pt\times\sum_{q=u,d,s}e_{q}^{2}\mathrm{Tr}\Big{(}\gamma_{\kappa}S_{q}(z,y)\gamma_{\sigma}S_{q}(y,z)$
$\displaystyle\hskip
56.9055pt-\big{\langle}\gamma_{\kappa}S_{q}(z,y)\gamma_{\sigma}S_{q}(y,z)\big{\rangle}_{\text{QCD}}\Big{)}\Big{\rangle}_{\text{QCD}},$

where $S_{q}$ denotes the quark propagator. As shown in Fig. 4, we calculate
the hadronic four-point function with two point-source propagators. The above
trick can make the evaluation more efficient [42, 47].

Figure 4: Two point-source propagators are used to calculate the hadronic
four-point function. The locations of the point-sources are indicated in the
diagrams as small circles.

We note that $\mathcal{H}^{\text{discon-no-perm}}$, the disconnected diagram
of the hadronic four-point function without permutation, still satisfies the
current conservation condition in the continuum limit (this is not true for
$\mathcal{H}^{\text{con-no-perm}}$) which leads to the following relation in
the infinite volume limit:

$\displaystyle\sum_{x_{\text{op}}}\mathcal{H}^{\text{discon-no-
perm}}_{k,\rho,\sigma,\lambda}(x_{\text{op}},x,y,z)=0.$ (25)

We are therefore allowed to alter the choice of $x_{\text{ref}}$ for the
disconnected diagram:

$\displaystyle x_{\text{ref-discon}}=x.$ (26)

This choice allows the summation over $x_{\text{op}}$ to be performed
independently of coordinates $y$ and $z$. This choice also has the benefit of
suppressing the contribution from the region where $|x-x_{\text{op}}|$ is
small and the quark loop is large before subtracting the vacuum expectation
value. Also, we note that the new choice $x_{\text{ref-discon}}$ is the same
as the initial choice in the long distance region where the contribution is
dominated by $\pi^{0}$-exchange and $|y-z|\ll\min(|x-y|,|x-z|)$. This property
guarantees that the connected and disconnected diagrams have exactly the same
QED weighting function in the long distance region, which leads to Eq. 42
before taking the continuum limit. This is the main reason for the choice of
$x_{\text{ref}}$ for the connected diagrams.

##  III Lattice details

The calculation presented here is performed on ensembles of gauge fields
generated by the RBC and UKQCD collaborations [52]. The main calculation is
carried out using the 48I ensemble, a physical mass ensemble generated with
2+1 flavors of Möbius domain wall fermions (MDWF). We also use a few other
ensembles to calculate various corrections and estimate systematic effects.
The relevant information about the 48I ensemble and other ensembles is listed
in Table 2. We always use the MDWF action. For most of the ensembles, the
quarks have their physical masses.

| 48I | 64I | 24D | 32D | 24DH
---|---|---|---|---|---
$m_{\pi}$ (MeV) | 139 | 135 | 142 | 142 | 341
$a^{-1}$ (GeV) | 1.730 | 2.359 | 1.015 | 1.015 | 1.015
$a$ (fm) | 0.114 | 0.084 | 0.194 | 0.194 | 0.194
$L$ (fm) | 5.47 | 5.38 | 4.67 | 6.22 | 4.67
$L_{s}$ | 24 | 12 | 24 | 24 | 24
# meas | 113 | 64 | 31 | 70 | 37
Table 1: 2+1 flavors of MDWF gauge field ensembles generated by the RBC-UKQCD
collaborations [52]. The labels indicate the lattice size in lattice units and
the QCD gauge action, where “I” stands for Iwasaki, and “D” stands for
Iwasaki+DSDR. DSDR stands for dislocation-suppressing-determinant-ratio, which
is an additional term in the action to soften explicit chiral symmetry
breaking effects, in particular for very coarse lattices [53]. The pion mass
$m_{\pi}$, lattice spacing $a$222After the analysis for this work was
completed, we updated the lattice spacing for the 24D and 32D ensembles. The
new lattice spacing is $a^{-1}=1.023~{}\mathrm{GeV}$ and makes no material
difference to the results presented here., spatial extent $L$, extra fifth
dimension size $L_{s}$, and the number of QCD configurations used are listed.

For our main calculation on the 48I ensemble, we use 113 configurations. On
each configuration, we randomly sample 2048 uniformly distributed distinct
points and calculate light quark point-source propagators for each of the
selected points.

Among the above 2048 points we randomly select 1024 points and also calculate
strange quark propagators. These 2048 + 1024 point-source propagators are
computed using the (deflated-) conjugate gradient algorithm with a sloppy
stopping condition. To further speed up the sloppy propagator calculation, we
use the zMöbius domain wall fermion formulation [54] with reduced fifth
dimension size $L_{s}=14$ to approximate the unitary MDWF action used in the
gauge ensemble generation ($L_{s}=24$). For deflation of the Dirac operator we
reuse the eigenvectors generated for the calculation of the hadronic vacuum
polarization contribution to the muon $g-2$ [55, 56] using the locally-
coherent Lanczos approach [57]. Then, we use a two-level all-mode-averaging
(AMA) method [58] to correct the bias caused by the inexact propagators. The
AMA method requires a portion of the sloppy propagators be computed again with
higher precision. In the first step a more precise version of the light and
strange quark propagators are calculated on a randomly selected set of 64
points among the 1024 points. In the second step among the 64 points, we
randomly select 16 points and calculate these propagators to full precision.
These propagators are used to calculate both the connected and the
disconnected diagrams in the bias correction. By using two steps, we are able
to reduce the error coming from the bias correction to a level well below the
statistical error for the sloppy part.

For the connected diagrams, we sample two-point pairs among all the possible
$2048\times(2048-1)/2+2048=2098176$ combinations (including the case where
$x=y$). For each of the two-point pairs, we can perform the contraction as
described in Eq. (II) in the previous section. As the number of all possible
pairs is enormous, the cost to perform all the contractions is not affordable.
Therefore, we only calculate the contraction for a subset of the available
two-point pairs. We sample the subset with the empirical probability $p(r)$,
which is a function of the distance between the two points, $r=|x-y|$,

$\displaystyle
p(r)=\left\\{\begin{array}[]{ll}\frac{N_{\text{psrc}}-1}{2(L^{3}T-1)}&\text{if
}r=0\\\ 1&\text{if }8\geq r>0\\\ \frac{1}{(r/8)^{8}}&\text{if }L\geq r>8\\\
0&\text{if }r>L\end{array}\right.,$ (31)

where $L=48$, $T=96$, and the numbers are in lattice units. On average, about
$57,000$ light quark point-source propagator pairs per configuration are
sampled to calculate the connected diagrams. The reason we sample long-
distance point-pairs less often is due to the fact that the connected hadronic
four point function and its statistical fluctuation decrease with distance for
long distances. For each pair, due to the sampling procedure, the probability
of that pair taking a particular relative coordinate is a function of $r$. The
total connected diagram contribution is equal to the expectation value of the
contraction for the point pairs multiplied by the inverse probabilities. We
refer to the inverse probability as the weight $w(r)$:

$\displaystyle w(r)=\left\\{\begin{array}[]{ll}w_{0}&\text{if }r=0\\\
w_{0}/p(r)&\text{if }r>0\end{array}\right.,$ (34) $\displaystyle
w_{0}=\sum_{x}p(|x|).$ (35)

To further reduce the contraction cost and the storage cost of saving these
propagators, we use field sparsening techniques [59, 60]. We randomly sample
$1/16$ points among all the points of the lattice and only perform
contractions on these randomly selected points. Since we only need propagator
values on these points, we only save these values to disk. To accommodate the
sparsening, we multiply a factor of $16^{2}$ for the summation over $z$ and
$x_{\text{op}}$ in Eq. (II), except when $z=x_{\text{op}}$, where we multiply
by $16$. In this way the sampling procedure does not introduce any systematic
effects on our central value.

For the disconnected diagrams, the two point-source locations are also
indicated in Fig. 4. Different from the connected diagrams, the statistical
fluctuations of the disconnected hadronic functions are not suppressed when
$r=|x-y|$ increases. Therefore, we use all possible combinations of point
pairs of $x$ and $y$ as long as $|x-y|\leq L$ to estimate the summation over
$x$ and $y$ in Eq. (II). Thanks to the choice of $x_{\text{ref}}=x_{\text{ref-
discon}}=x$ in Eq. (26), the summation over $x_{\text{op}}$ can be performed
for each point-source propagator with source location at $x$, the result can
be used for all possible values of $y$ locations. However, the summation of
$z$ depends on the location of $x$ and $y$ due to the QED weighting function.
Therefore, we need to perform the summation over $z$ for each pair of points.

To make the contractions of all point pairs affordable, we aggressively
sparsen when summing over $z$. Fortunately, the vertex $z$ and vertex $y$ are
connected by two quark propagators, and the four point function is
exponentially suppressed when the distance between $z$ and $y$ increases.

Naturally, we can sample the $z$ locations based on the distance between $z$
and $y$, similar to the sampling of $x$ and $y$ combinations in the
calculation of the connected diagrams. However, for this disconnected diagram,
we discovered a much more efficient adaptive sampling scheme as follows. For
each point-source location $y$, we calculate the following square norm
$n(z,y)$ of the quark loop for all $z$:

$\displaystyle n(z,y)$
$\displaystyle=\sum_{\kappa,\sigma}\Big{|}\mathrm{Tr}\Big{(}\gamma_{\kappa}S_{q}(z,y)\gamma_{\sigma}S_{q}(y,z)$
$\displaystyle\hskip
28.45274pt-\big{\langle}\gamma_{\kappa}S_{q}(z,y)\gamma_{\sigma}S_{q}(y,z)\big{\rangle}_{\text{QCD}}\Big{)}\Big{|}^{2}$

Unlike the full contraction, the above norm is independent of the location of
$x$, and the location of $z$ is sampled with the following (empirical)
probability

$\displaystyle p_{y}(z)=\left\\{\begin{array}[]{ll}1&\text{if }n(z,y)\geq
t_{0}^{2}\text{ and }|z-y|\leq L\\\ \sqrt{n(z,y)}/t_{0}&\text{if
}n(z,y)<t_{0}^{2}\text{ and }|z-y|\leq L\\\ 0&\text{if
}|z-y|>L\end{array}\right.,$ (39)

where $t_{0}$ is the threshold parameter to control the overall sample
frequency. With this probability distribution, we sample more often points
where the (norm of) the quark loop is large which is much more efficient than
basing the probability on the distance $|z-y|$. We choose $t_{0}=10^{-3}$ in
this work for the 48I ensemble.

For the light quark loop, there are $9,156,431$ points within the $|z-y|\leq
L$ range. Among these, we sample about $17,000$ points for $z$ on average. For
each of the sampled $z$, we then loop over all possible $x$ locations with
$|x-y|\leq L$ and perform the full contraction. Due to the sampling procedure,
we multiply the result of the contraction by the inverse of the sample
probability $1/p_{y}(z)$ to obtain an unbiased final result. While the sampled
$z$ points amount to less than $0.2\%$ of the possible $z$ points in the
allowed range, we expect the final statistical precision to be almost
unaffected by the adaptive sampling procedure. That is, we expect the final
statistical error would be almost the same even if we had calculated the
contraction using all the $z$ locations. This expectation is based on a quick
test run with much lower threshold $t_{0}=5\times 10^{-5}$, which corresponds
to about a factor of $20$ smaller number of sampled points; we find the final
statistical error is roughly the same.

In contrast to the connected diagram, we do not sparsen the propagators. To
save disk space, we temporarily save the following intermediate contraction
for each point-source propagator:

$\displaystyle\mathrm{Tr}\Big{(}\gamma_{\kappa}S_{q}(z,y)\gamma_{\sigma}S_{q}(y,z)\Big{)}.$
(40)

For each propagator on each lattice site, the above contraction corresponds to
$4\times 4$ complex numbers while the full propagator takes $12\times 12$
complex numbers.

##  IV Results

In this section we display several figures of the summand (integrand) or its
partial-sum as a function of $R_{\text{max}}$:

$\displaystyle R_{\text{max}}=\text{max}(|x-y|,|y-z|,|x-z|),$ (41)

where $x$, $y$, and $z$ represent the positions of the three internal quark-
photon vertices. Partial-sums are obtained by summing the corresponding
summand from $0$ up to the specified value of $R_{\text{max}}$ (inclusive).
The rightmost value in a partial-sum plot corresponds to the total
contribution.

Since the summand decreases exponentially at large $R_{\text{max}}$, we expect
the partial-sum to approach a plateau for large enough $R_{\text{max}}$. We
usually use the result of the partial-sum at $4~{}\mathrm{fm}$, which includes
all contributions from the region $R_{\text{max}}<4~{}\mathrm{fm}$. In the
lattice calculation, we calculate the partial sum with respect to
$R_{\text{max}}$ and save the partial sum for all integer values (in lattice
units) of the upper limit of $R_{\text{max}}$. To compare with results
calculated using ensembles with different lattice spacings, we can linearly
interpolate the partial sum for any upper limits of $R_{\text{max}}$. In this
work, we choose to interpolate at integer multiples
$\times\,0.1~{}\mathrm{fm}$. The summand is obtained at half-integer multiples
$\times\,0.1~{}\mathrm{fm}$ by taking the difference of the linearly
interpolated partial sum. Therefore, the data points of the plots shown in
this work always have spacing $0.1~{}\mathrm{fm}$, and the data points of the
integrand plots always have its x-axis value equal to half-integer times
$0.1~{}\mathrm{fm}$. In this work, we present our results of anomalous
magnetic moment in the unit of $10^{-10}$.

###  IV.1 Light quark contribution

Figure 5: Light quark contributions computed on the 48I ensemble from the
connected diagrams, the disconnected diagrams, and the total. The upper plot
shows the corresponding summands and the lower plot shows the partial sum.

Figure 6: Light quark contributions from a special combination of the
connected and disconnected diagrams to cancel the long-distance $\pi^{0}$
exchange contribution (see Eq. 43). The upper plot shows the summand and the
lower plot shows the partial sum. The combined summand vanishes much faster
than the light quark diagrams by themselves. Solid curves represent a fit of
the data to Eq. (44). The fit starts at $0.5~{}\mathrm{fm}$. The dashed lines
denote the statistical uncertainty of the fit.

The contributions from the connected diagrams, the disconnected diagrams, and
the sum of the two contributions are shown in Fig. 5. As can be seen from the
figure, the statistical error grows as $R_{\text{max}}$ increases. However,
for both the connected and the disconnected diagrams, there are non-negligible
contributions that come from regions where $R_{\text{max}}$ is as large as 3
or 4 fm. The statistical error is quite significant. Making the situation
worse, the contributions from the long-distance region of the connected and
the disconnected diagrams are opposite in sign, while the statistical error of
the connected and the disconnected diagrams are largely independent. As a
result, the relative uncertainty on the sum is much larger than that for the
individual contributions. Results are given in Tab. 2 as “48I light con
$R_{\text{max}}<4\mathrm{fm}$”, “48I light discon
$R_{\text{max}}<4\mathrm{fm}$”, and “48I light $R_{\text{max}}<4\mathrm{fm}$”.

To improve the situation, note that the large contribution and the
cancellation between the connected and disconnected diagrams at long distance
are due to $\pi^{0}$ exchange and are well-understood theoretically [61, 62,
48]. It has been shown that, at large $R_{\text{max}}$, the ratio of the
disconnected and the connected hadronic four-point function is $-25/34$. In
our present computational setup, we use the same infinite volume QCD weighting
function for both the connected and disconnected diagrams, and use the same
variable $R_{\text{max}}$ to study the partial-sum of the connected and
disconnected diagrams.333One can use different setup for the connected and
disconnected diagrams, in which case this relation will not hold exactly. For
example, Ref. [48] uses different setups. Therefore the same ratio applies to
the contribution to $a_{\mu}$. Formally, we have:

$\displaystyle\lim_{R\to\infty}\frac{a_{\mu}^{\text{discon}}(R_{\text{max}}>R)}{a_{\mu}^{\text{con}}(R_{\text{max}}>R)}=-\frac{1}{2}\cdot\frac{(e_{u}^{2}+e_{d}^{2})^{2}}{e_{u}^{4}+e_{d}^{4}}=-\frac{25}{34}.$
(42)

This ratio is exact in the large $R_{\text{max}}$ limit and not affected by
lattice artifacts or finite volume effects. Therefore, we can construct the
following combination:

$\displaystyle a_{\mu}^{\text{no-
pion}}=a_{\mu}^{\text{discon}}+\frac{25}{34}a_{\mu}^{\text{con}},$ (43)

where the $\pi^{0}$ exchange contribution to $a_{\mu}^{\text{no-pion}}$
vanishes in the long distance. We plot the summand and partial-sum of
$a_{\mu}^{\text{no-pion}}$ in Fig. 6. Indeed, the partial sum of
$a_{\mu}^{\text{no-pion}}$ reaches the plateau much earlier than the connected
or disconnected diagrams. This trick of combining the connected and
disconnected with appropriate factors to obtain a faster plateau was employed
in Ref. [64]. For $a_{\mu}^{\text{no-pion}}$, we will use $2.0~{}\mathrm{fm}$
and $2.5~{}\mathrm{fm}$ as the upper limit of $R_{\text{max}}$. The results
are shown in Tab. 2 as “48I light no-pion $R_{\text{max}}<2.5\mathrm{fm}$” or
“48I light no-pion $R_{\text{max}}<2\mathrm{fm}$”.

In the upper panel of Fig. 6, we fit the summand to the following empirical
form:

$\displaystyle
f(R_{\text{max}})=A\frac{R_{\text{max}}^{6}}{R_{\text{max}}^{3}+C^{3}}e^{-BR_{\text{max}}}$
(44)

The fit range is from 0.5 fm to 4 fm. We use the result of the fit to estimate
the long-distance contribution to $a_{\mu}^{\text{no-pion}}$. Since this is a
completely empirical fit, we will assign a $100\%$ systematic uncertainty to
it. The “no-pion" results are also collected in Tab. 2.

Contribution name | $a_{\mu}$ | $\times 10^{10}$
---|---|---
48I light con $R_{\text{max}}<2\mathrm{fm}$ | $7.24$ | $(0.15)_{\text{stat}}$
48I light con $R_{\text{max}}<2.5\mathrm{fm}$ | $10.64$ | $(0.31)_{\text{stat}}$
48I light con $R_{\text{max}}<4\mathrm{fm}$ | $18.61$ | $(1.22)_{\text{stat}}$
48I light con $R_{\text{max}}>4\mathrm{fm}$ | $7.56$ | $(0.42)_{\text{stat}}(1.06)_{\text{syst}}~{}[1.14]$
48I light con FV-corr | $-1.79$ | $\phantom{(0.00)_{\text{stat}}}(0.42)_{\text{syst}}$
48I light con $m_{\pi}$-corr | $1.32$ | $(0.27)_{\text{stat}}(0.66)_{\text{syst}}~{}[0.72]$
48I light con $a^{2}$-corr | $0.00$ | $\phantom{(0.00)_{\text{stat}}}(1.49)_{\text{syst}}$
light con | $25.70$ | $(1.33)_{\text{stat}}(1.99)_{\text{syst}}~{}[2.39]$
48I light no-pion $R_{\text{max}}<2.5\mathrm{fm}$ | $5.09$ | $(0.76)_{\text{stat}}$
48I light no-pion $R_{\text{max}}<2\mathrm{fm}$ | $4.65$ | $(0.42)_{\text{stat}}$
48I light no-pion $R_{\text{max}}>2.5\mathrm{fm}$ | $0.31$ | $(0.22)_{\text{stat}}(0.31)_{\text{syst}}~{}[0.38]$
48I light no-pion $R_{\text{max}}>2\mathrm{fm}$ | $0.88$ | $(0.46)_{\text{stat}}(0.88)_{\text{syst}}~{}[0.99]$
48I light discon $R_{\text{max}}<2\mathrm{fm}$ | $-0.67$ | $(0.41)_{\text{stat}}$
48I light discon $R_{\text{max}}<2.5\mathrm{fm}$ | $-2.73$ | $(0.74)_{\text{stat}}$
48I light discon $R_{\text{max}}<4\mathrm{fm}$ | $-7.49$ | $(1.82)_{\text{stat}}$
48I light discon $R_{\text{max}}<4\mathrm{fm}$ hybrid-2.5fm | $-8.28$ | $(1.31)_{\text{stat}}(0.31)_{\text{syst}}~{}[1.35]$
48I light discon $R_{\text{max}}<4\mathrm{fm}$ hybrid-2fm | $-8.15$ | $(1.24)_{\text{stat}}(0.88)_{\text{syst}}~{}[1.51]$
48I light discon $R_{\text{max}}>4\mathrm{fm}$ | $-5.56$ | $(0.31)_{\text{stat}}(0.78)_{\text{syst}}~{}[0.84]$
48I light discon FV-corr | $1.31$ | $\phantom{(0.00)_{\text{stat}}}(0.31)_{\text{syst}}$
48I light discon $m_{\pi}$-corr | $-0.98$ | $(0.20)_{\text{stat}}(0.49)_{\text{syst}}~{}[0.53]$
48I light discon $a^{2}$-corr | $0.00$ | $\phantom{(0.00)_{\text{stat}}}(0.66)_{\text{syst}}$
light discon | $-12.71$ | $(1.87)_{\text{stat}}(1.17)_{\text{syst}}~{}[2.20]$
light discon hybrid-2.5fm | $-13.50$ | $(1.36)_{\text{stat}}(1.21)_{\text{syst}}~{}[1.82]$
light discon hybrid-2fm | $-13.37$ | $(1.29)_{\text{stat}}(1.46)_{\text{syst}}~{}[1.95]$
48I light $R_{\text{max}}<2\mathrm{fm}$ | $6.57$ | $(0.43)_{\text{stat}}$
48I light $R_{\text{max}}<2.5\mathrm{fm}$ | $7.90$ | $(0.78)_{\text{stat}}$
48I light $R_{\text{max}}<4\mathrm{fm}$ | $11.11$ | $(2.11)_{\text{stat}}$
48I light $R_{\text{max}}<4\mathrm{fm}$ hybrid-2.5fm | $10.32$ | $(0.99)_{\text{stat}}(0.31)_{\text{syst}}~{}[1.04]$
48I light $R_{\text{max}}<4\mathrm{fm}$ hybrid-2fm | $10.46$ | $(0.89)_{\text{stat}}(0.88)_{\text{syst}}~{}[1.25]$
48I light $R_{\text{max}}>4\mathrm{fm}$ | $2.00$ | $(0.11)_{\text{stat}}(0.28)_{\text{syst}}~{}[0.30]$
48I light FV-corr | $-0.47$ | $\phantom{(0.00)_{\text{stat}}}(0.11)_{\text{syst}}$
48I light $m_{\pi}$-corr | $0.35$ | $(0.07)_{\text{stat}}(0.17)_{\text{syst}}~{}[0.19]$
48I light $a^{2}$-corr | $0.00$ | $\phantom{(0.00)_{\text{stat}}}(0.83)_{\text{syst}}$
light total | $12.99$ | $(2.11)_{\text{stat}}(0.90)_{\text{syst}}~{}[2.29]$
light total hybrid-2.5fm | $12.20$ | $(1.01)_{\text{stat}}(0.95)_{\text{syst}}~{}[1.38]$
light total hybrid-2fm | $12.33$ | $(0.90)_{\text{stat}}(1.25)_{\text{syst}}~{}[1.55]$
Table 2: Light quark contributions computed on the 48I ensemble. Values are
different contributions to $a_{\mu}\times 10^{10}$ where
$a_{\mu}=(g_{\mu}-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.

With the definition of $a_{\mu}^{\text{no-pion}}$, we obtain:

$\displaystyle a_{\mu}^{\text{discon}}$ $\displaystyle=a_{\mu}^{\text{no-
pion}}-\frac{25}{34}a_{\mu}^{\text{con}},$ (45) $\displaystyle
a_{\mu}^{\text{total}}$ $\displaystyle=a_{\mu}^{\text{no-
pion}}+\frac{9}{34}a_{\mu}^{\text{con}}.$ (46)

The advantage comes in using the early plateau value of $a_{\mu}^{\text{no-
pion}}$ and combining it with $a_{\mu}^{\text{con}}$, which plateaus at much
larger $R_{\text{max}}$. This is a “hybrid” approach to calculate
$a_{\mu}^{\text{discon}}$ and $a_{\mu}^{\text{total}}$, which has a much
smaller statistical error than the direct combination. For
$a_{\mu}^{\text{con}}$, we will use $4~{}\mathrm{fm}$ as the upper limit of
$R_{\text{max}}$. The results are shown in Tab. 2 as entries that start with
“48I light discon $R_{\text{max}}<4\mathrm{fm}$ hybrid-” and “48I light
$R_{\text{max}}<4\mathrm{fm}$ hybrid-”.

The fitting function form in Eq. (44) is inspired by the function in Eq. (20)
of Ref. [48]. The functional form used here differs due to the meanings of
$R_{\text{max}}$ and the variable “$|y|$” in Ref. [48]. Also, note that the
subtraction scheme of the QED weighting function is also different [46, 47,
48, 65]. We tried using this fit function to fit other contributions. The
results are shown in Appendix C.

The remaining finite volume, pion mass, and non-zero lattice spacing
corrections will be studied in Sections IV.4, IV.5, and IV.6.

###  IV.2 Strange quark contribution

Figure 7: Strange quark contributions computed on the 48I ensemble from the
connected diagrams, the disconnected diagrams, and the total. The upper plot
shows the corresponding summands and the lower plot shows the partial sum.

The contributions from the strange quark connected diagrams, the disconnected
diagrams, and the sum of the two contributions are plotted in Fig. 7. Note the
strange quark disconnected diagrams include diagrams where one or both loops
are strange quark loops, while the light quark disconnected diagrams discussed
in the previous section only contain light quarks. As seen in the figure, the
contribution from the strange quark-connected diagrams is very precise and
very small. It can also be clearly seen that the strange quark contribution
vanishes much faster at long-distance compared to the light quark
contribution.

The disconnected diagrams, on the other hand, are much noisier. However, we
still expect the signal vanishes faster than the light quark contribution due
to the absence of the $\pi^{0}$ exchange contribution. Therefore, we can treat
the strange quark disconnected diagrams similar to $a_{\mu}^{\text{no-pion}}$.
We listed the results with $4.0~{}\mathrm{fm}$, $2.5~{}\mathrm{fm}$, and
$2~{}\mathrm{fm}$ as the upper limit of $R_{\text{max}}$. These results are
shown in Tab. 3 as entries that start with “48I strange discon
$R_{\text{max}}<$”. We estimate the systematic error caused by truncating the
strange quark disconnected diagram integration using the size of the long-
distance contribution from the strange quark connected diagrams, which is also
shown in Tab. 3 as entries that start with “48I strange con
$R_{\text{max}}>$”. The size of the long-distance connected contribution is
used as the estimation for the systematic uncertainty for the corresponding
disconnected contributions.

The remaining finite volume, pion mass, and non-zero lattice spacing
corrections will be studied in Section. IV.4, IV.5, and IV.6.

Contribution name | $a_{\mu}$ | $\times 10^{10}$
---|---|---
48I strange con $R_{\text{max}}<4\mathrm{fm}$ | $0.327$ | $(0.002)_{\text{stat}}$
48I strange con $a^{2}$-corr | $0.026$ | $(0.008)_{\text{stat}}$
strange con | $0.353$ | $(0.007)_{\text{stat}}$
48I strange con $R_{\text{max}}>2.5\mathrm{fm}$ | $0.006$ | $(0.000)_{\text{stat}}$
48I strange con $R_{\text{max}}>2\mathrm{fm}$ | $0.024$ | $(0.000)_{\text{stat}}$
48I strange discon $R_{\text{max}}<4\mathrm{fm}$ | $-0.380$ | $(0.607)_{\text{stat}}$
48I strange discon $R_{\text{max}}<2.5\mathrm{fm}$ | $-0.357$ | $(0.223)_{\text{stat}}$
48I strange discon $R_{\text{max}}<2\mathrm{fm}$ | $-0.280$ | $(0.121)_{\text{stat}}$
48I strange discon $a^{2}$-corr | $0.000$ | $\phantom{(0.000)_{\text{stat}}}(0.029)_{\text{syst}}$
strange discon | $-0.380$ | $(0.607)_{\text{stat}}(0.029)_{\text{syst}}~{}[0.608]$
strange discon hybrid-2.5fm | $-0.357$ | $(0.223)_{\text{stat}}(0.029)_{\text{syst}}~{}[0.225]$
strange discon hybrid-2fm | $-0.280$ | $(0.121)_{\text{stat}}(0.037)_{\text{syst}}~{}[0.126]$
48I strange $R_{\text{max}}<4\mathrm{fm}$ | $-0.053$ | $(0.607)_{\text{stat}}$
48I strange $R_{\text{max}}<4\mathrm{fm}$ hybrid-2.5fm | $-0.030$ | $(0.222)_{\text{stat}}(0.006)_{\text{syst}}~{}[0.223]$
48I strange $R_{\text{max}}<4\mathrm{fm}$ hybrid-2fm | $0.047$ | $(0.121)_{\text{stat}}(0.024)_{\text{syst}}~{}[0.123]$
48I strange $a^{2}$-corr | $0.026$ | $(0.008)_{\text{stat}}(0.029)_{\text{syst}}~{}[0.030]$
strange total | $-0.027$ | $(0.607)_{\text{stat}}(0.029)_{\text{syst}}~{}[0.608]$
strange total hybrid-2.5fm | $-0.004$ | $(0.223)_{\text{stat}}(0.029)_{\text{syst}}~{}[0.225]$
strange total hybrid-2fm | $0.073$ | $(0.121)_{\text{stat}}(0.037)_{\text{syst}}~{}[0.127]$
Table 3: Strange quark contributions computed on the 48I ensemble. Values are
different contributions to $a_{\mu}\times 10^{10}$ where
$a_{\mu}=(g_{\mu}-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.

###  IV.3 Long-distance $\pi^{0}$-exchange contribution

Figure 8: The long-distance HLbL contribution to the muon $g-2$ associated
with $\pi^{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-volume, continuum expression.

The dominant source of the HLbL contribution is expected to be the so-called
“$\pi^{0}$-pole” contribution [5]. Due to the small mass of the $\pi^{0}$,
this contribution can be non-negligible for large $R_{\text{max}}$ in a
Euclidean space-time lattice calculation, even where long-distance
contributions and finite volume effects are exponentially suppressed since it
is only suppressed as $e^{-m_{\pi}R_{\text{max}}}$.

We separately calculate the long-distance $\pi^{0}$-exchange contribution as
illustrated in Fig. 8. We use the lattice calculation of the $\pi^{0}$
transition form factor:

$\displaystyle\mathcal{F}_{\mu,\nu}(x,p)$ $\displaystyle=\langle
0|TJ_{\mu}(x)J_{\nu}(0)|\pi^{0}(\vec{p})\rangle$ (47)
$\displaystyle=\epsilon_{\mu,\nu,\rho,\sigma}x_{\rho}p_{\sigma}\mathcal{F}(x^{2},p\cdot
x),$ (48)

where $p$ is the Euclidean four-momentum of the $\pi^{0}$ state. We use
$p=(im_{\pi},\vec{0})=im_{\pi}\hat{t}$ in the numerical lattice calculation.
So we can obtain the numerical values of

$\displaystyle\mathcal{F}_{\mu,\nu}\big{(}x,im_{\pi}\hat{t}\big{)}$ (49)

for all possible Euclidean space-time locations $x$ directly from lattice QCD
calculations, without assuming any particular function form for the transition
form factor. This form factor can then be used to construct approximations for
infinite volume hadronic correlation functions with large separation. A
detailed derivation is given in Appendix B. Here, we briefly describe the
idea. First, we introduce a properly normalized $\pi^{0}$ interpolating field,
e.g.

$\displaystyle\pi^{0}(x)=Z_{\pi}^{1/2}\frac{i}{\sqrt{2}}\big{(}\bar{u}(x)\gamma_{5}u(x)-\bar{d}(x)\gamma_{5}d(x)\big{)}\,.$
(50)

The normalization constant $Z_{\pi}$ is determined by the following
requirement:

$\displaystyle\langle 0|\pi^{0}(x)|\pi^{0}(\vec{p})\rangle=e^{ip\cdot x}.$
(51)

We can then rewrite Eq. (47) in a slightly different form:

$\displaystyle\langle 0|TJ_{\mu}(x)J_{\nu}(y)|\pi^{0}(\vec{p})\rangle$
$\displaystyle\qquad=\mathcal{F}_{\mu,\nu}(x-y,p)\langle
0|\pi^{0}(y)|\pi^{0}(\vec{p})\rangle\,.$ (52)

The above relation suggest that the timer-ordered product
$TJ_{\mu}(x)J_{\nu}(y)$, when used in between the vacuum state and a single
$\pi^{0}$ state, can be viewed as a properly weighted $\pi^{0}$ interpolating
field: $\mathcal{F}_{\mu,\nu}(x-y,p)\pi^{0}(y)$. This property can be used to
calculate the following three-point function:

$\displaystyle\langle TJ_{\mu}(x)J_{\nu}(y)\pi^{0}(z)\rangle$ (53)

For very large $m_{\pi}|y-z|$ and $y-z$ along the positive time direction,
i.e. $y-z=(|y-z|,\vec{0})$, the intermediate states between
$TJ_{\mu}(x)J_{\nu}(y)$ and $\pi^{0}(z)$ will be mostly $\pi^{0}$ states with
$\vec{p}\approx\vec{0}$. We can employ Eq. (IV.3) with
$p\approx(im_{\pi},\vec{0})$ and obtain:

$\displaystyle\langle TJ_{\mu}(x)J_{\nu}(y)\pi^{0}(z)\rangle$
$\displaystyle\approx\mathcal{F}_{\mu,\nu}\Big{(}x-y,im_{\pi}\frac{y-z}{|y-z|}\Big{)}\langle
T\pi^{0}(y)\pi^{0}(z)\rangle.$ (54)

Note that we have expressed the above relation in an $O(4)$ rotationally
covariant form, so the equation is also valid for $y-z$ in directions other
than the time direction. Then, we can introduce the infinite-volume free
scalar propagator $D_{\pi^{0}}(z)$, with physical pion mass. For sufficiently
large $|x-y|$, we have:

$\displaystyle\langle T\pi^{0}(x)\pi^{0}(y)\rangle=D_{\pi^{0}}(x-y).$ (55)

Now, the relation between Eq. (47) and Eq. (IV.3) is established. We can
approximate the infinite volume hadronic four-point function illustrated in
Fig. 8 in the large $|x-y|$ region in a similar approach:

$\displaystyle\langle
TJ_{\mu^{\prime}}(x^{\prime})J_{\mu}(x)J_{\nu^{\prime}}(y^{\prime})J_{\nu}(y)\rangle\approx
D_{\pi^{0}}(x-y)$
$\displaystyle\times\mathcal{F}_{\mu^{\prime},\mu}\Big{(}x^{\prime}-x,im_{\pi}\frac{x-y}{|x-y|}\Big{)}$
$\displaystyle\times\mathcal{F}_{\nu^{\prime},\nu}\Big{(}y^{\prime}-y,im_{\pi}\frac{y-x}{|y-x|}\Big{)}.$
(56)

In the above approximation, $|x-x^{\prime}|$ and $|y-y^{\prime}|$ will not be
too large in order to create localized pion states, so the form factor
$\mathcal{F}$ can be computed using a reasonably sized lattice. As described
in Eq. (47) and Eq. (49), we only directly calculate the transition form
factor $\mathcal{F}$ for a zero-momentum pion state. To obtain the form factor
needed in Eq. (IV.3) above, we need to perform $O(4)$ rotations:

$\displaystyle\mathcal{F}_{\mu,\nu}\Big{(}\widetilde{x},im_{\pi}\hat{n}\Big{)}=\Lambda_{\mu,\mu^{\prime}}\Lambda_{\nu,\nu^{\prime}}\mathcal{F}_{\mu^{\prime},\nu^{\prime}}\Big{(}\widetilde{x}^{\prime},im_{\pi}\hat{t}\Big{)},$
(57)

where $\tilde{x}$ can be either $x^{\prime}-x$ or $y^{\prime}-y$ in Eq. (IV.3)
and $\hat{n}$ can be $\pm(x-y)/|x-y|$. The $O(4)$ rotation matrix
$\Lambda=\Lambda(\hat{n})$ and Euclidean space-time coordinate
$\widetilde{x}^{\prime}$ satisfy:

$\displaystyle\delta_{\mu,\nu}$
$\displaystyle=\Lambda_{\mu,\mu^{\prime}}\Lambda_{\nu,\nu^{\prime}}\delta_{\mu^{\prime},\nu^{\prime}},$
(58) $\displaystyle\hat{n}_{\mu}$
$\displaystyle=\Lambda_{\mu,\mu^{\prime}}\hat{t}_{\mu^{\prime}},$ (59)
$\displaystyle\widetilde{x}_{\mu}$
$\displaystyle=\Lambda_{\mu,\mu^{\prime}}\widetilde{x}^{\prime}_{\mu^{\prime}}.$
(60)

We would like to emphasize again that we do not assume any particular form of
the $\pi^{0}$ transition form factor and perform fits. The inputs to Eq.
(IV.3) can be directly obtained from the $O(4)$ rotation in Eq. (57) and the
lattice QCD calculation of the position-space matrix elements
$\mathcal{F}_{\mu,\nu}(x,im_{\pi}\hat{t})$ for a zero-momentum $\pi^{0}$ state
defined in Eq. (47) and Eq. (49). The only approximation is made in Eq.
(IV.3), which is based on large $m_{\pi}|x-y|$.

Using the 48I ensemble to calculate $\mathcal{F}$ and making this
approximation, we obtained the long-distance part corresponding to
$R_{\text{max}}>4~{}\mathrm{fm}$:

$\displaystyle a_{\mu}(R_{\text{max}}>4~{}\mathrm{fm})\times 10^{10}$ (61)
$\displaystyle\approx
a_{\mu}^{\pi^{0}\text{-exch}}(R_{\text{max}}>4~{}\mathrm{fm})\times
10^{10}=2.00(11)_{\text{stat}}(28)_{\text{syst}},$

where we use the same subtracted infinite-volume QED weighting function as in
the direct calculation. We estimate the relative systematic uncertainty from
the approximation employed in Eq. (IV.3) due to not including the charged loop
contribution as:

$\displaystyle\frac{e^{-2m_{\pi}(4\,\mathrm{fm})}}{e^{-m_{\pi}(4\,\mathrm{fm})}}\approx
6\%$ (62)

and the error from approximating the $\pi^{0}$ propagation direction to be
along $x-y$ to be:

$\displaystyle\frac{0.5\,\mathrm{fm}}{4\,\mathrm{fm}}\approx 13\%$ (63)

where we assume the typical separation for $|x-x^{\prime}|$ and
$|y-y^{\prime}|$ to be $0.5~{}\mathrm{fm}$ and $|x-y|$ to be
$4~{}\mathrm{fm}$. Combining these two estimates in quadrature, we assigned
14% total systematic uncertainty for the long distance part
($R_{\text{max}}>4~{}\mathrm{fm}$) contribution. Using the ratio between the
connected and disconnected contributions as described in Eq. 42, we can also
properly assign this contribution to the connected and disconnected diagrams.
The results are shown in Tab. 2 with the labels “$\dots
R_{\text{max}}>4\mathrm{fm}$”.

###  IV.4 Finite volume corrections

The long-distance part of the HLbL contribution usually suffers the most
significant finite volume effects. However, the long-distance $\pi^{0}$
exchange contribution calculated in Section IV.3 is performed in the infinite
volume and is free of finite volume effects. Therefore, we only need to
correct for the finite volume effects for the relatively short-distance
region, where $R_{\text{max}}<4~{}\mathrm{fm}$. These finite volume effects
are expected to be quite small due to the constraint
$R_{\text{max}}<4~{}\mathrm{fm}$. The finite volume effects can be estimated
with the “$\pi^{0}$-pole” contribution as defined in Ref. [20]. Similar to
Ref. [66, 21], we define the $\pi^{0}$ transition form factor in Euclidean
space-time as:

$\displaystyle\mathcal{F}_{\mu,\nu}(x,p)=\langle
0|TJ_{\mu}(x)J_{\nu}(0)|\pi^{0}(\vec{p})\rangle$ $\displaystyle=$
$\displaystyle\int\frac{d^{4}q_{1}}{(2\pi)^{4}}e^{iq_{1}\cdot
x}\frac{-i}{4\pi^{2}F_{\pi}}\epsilon_{\mu,\nu,\rho,\sigma}{q_{1}}_{\rho}{q_{2}}_{\sigma}F_{\pi^{0}\gamma\gamma}(q_{1}^{2},q_{2}^{2})$
(64)

where $p=q_{1}+q_{2}$ and $p^{2}=-m_{\pi}^{2}$. We can convert the form factor
into coordinate space via the following:

$\displaystyle\int d^{4}u\,e^{ip\cdot
u}\widetilde{\mathcal{F}}_{\mu,\nu}(u,x,y)=\mathcal{F}_{\mu,\nu}(x-y,p).$ (65)

The solution to the above condition is not unique. Based on the momentum space
form factor $F_{\pi^{0}\gamma\gamma}(q_{1}^{2},q_{2}^{2})$ in Eq. (IV.4), we
obtain

$\displaystyle\widetilde{\mathcal{F}}_{\mu,\nu}(u,x,y)$
$\displaystyle=\frac{i}{4\pi^{2}F_{\pi}}\epsilon_{\mu,\nu,\rho,\sigma}\partial^{x}_{\rho}\partial^{y}_{\sigma}$
(66)
$\displaystyle\times\int\frac{d^{4}q_{1}}{(2\pi)^{4}}e^{iq_{1}\cdot(x-u)}\int\frac{d^{4}q_{2}}{(2\pi)^{4}}e^{iq_{2}\cdot(y-u)}F_{\pi^{0}\gamma\gamma}(q_{1}^{2},q_{2}^{2}).$

Note this definition of the position space form factor is different from the
inverse Fourier transformation of the Euclidean space hadronic matrix element
$\mathcal{F}_{\mu,\nu}(x,p)$, since the on-shell condition
$p^{2}=-m_{\pi}^{2}$ for the physical $\pi^{0}$ state cannot be satisfied
always in the inverse Fourier transformation. With the form factor defined
above, we can construct the $\pi^{0}$-pole contribution to the hadronic four-
point function:

$\displaystyle\langle
TJ_{\mu}(x)J_{\mu^{\prime}}(x^{\prime})J_{\nu}(y)J_{\nu^{\prime}}(y^{\prime})\rangle$
(67) $\displaystyle\approx\int d^{4}u\int d^{4}v\,D_{\pi^{0}}(u-v)$
$\displaystyle\times\Big{(}\widetilde{\mathcal{F}}_{\mu,\mu^{\prime}}(u,x,x^{\prime})\widetilde{\mathcal{F}}_{\nu,\nu^{\prime}}(v,y,y^{\prime})$
$\displaystyle\hskip
28.45274pt+\widetilde{\mathcal{F}}_{\mu,\nu}(u,x,y)\widetilde{\mathcal{F}}_{\mu^{\prime},\nu^{\prime}}(v,x^{\prime},y^{\prime})$
$\displaystyle\hskip
28.45274pt+\widetilde{\mathcal{F}}_{\mu,\nu^{\prime}}(u,x,y^{\prime})\widetilde{\mathcal{F}}_{\mu^{\prime},\nu}(v,x^{\prime},y)\Big{)}$

where $D_{\pi^{0}}(x-y)$ is the infinite volume free scalar propagator with
physical pion mass. In this calculation, we use the Lowest Meson Dominance
(LMD) model [38, 67, 68] for the $\pi^{0}$ transition form factors.

$\displaystyle F_{\pi^{0}\gamma\gamma}(q_{1}^{2},q_{2}^{2})\approx
F^{\text{LMD}}_{\pi^{0}\gamma\gamma}(q_{1}^{2},q_{2}^{2})$ $\displaystyle=$
$\displaystyle
F^{\text{VMD}}_{\pi^{0}\gamma\gamma}(q_{1}^{2},q_{2}^{2})+\frac{8\pi^{2}F_{\pi}^{2}}{3m_{V}^{2}}$
(68) $\displaystyle\hskip
5.69046pt\times\Big{(}F^{\text{TE}}_{\pi^{0}\gamma\gamma}(q_{1}^{2},q_{2}^{2})-F^{\text{VMD}}_{\pi^{0}\gamma\gamma}(q_{1}^{2},q_{2}^{2})\Big{)},$

where

$\displaystyle
F^{\text{VMD}}_{\pi^{0}\gamma\gamma}(q_{1}^{2},q_{2}^{2})=\frac{m_{V}^{2}}{q_{1}^{2}+m_{V}^{2}}\frac{m_{V}^{2}}{q_{2}^{2}+m_{V}^{2}}$
(69) $\displaystyle
F^{\text{TE}}_{\pi^{0}\gamma\gamma}(q_{1}^{2},q_{2}^{2})=\frac{m_{V}^{2}/2}{q_{1}^{2}+m_{V}^{2}}+\frac{m_{V}^{2}/2}{q_{2}^{2}+m_{V}^{2}}.$
(70)

The parameters used in the calculation are $m_{\pi}=134.9766~{}\mathrm{MeV}$,
$F_{\pi}=92~{}\mathrm{MeV}$, $m_{V}=770~{}\mathrm{MeV}$. We discretize the
model and calculate it with lattices of different sizes using the same
subtracted infinite volume QED weighting as the direct calculation. The
results are given in Tab. 4. Note the physical size of the “Match” and
“Coarse” lattices are the same as the 48I ensemble, which we used to perform
our main lattice QCD calculation. The spatial size is $L=5.5~{}\mathrm{fm}$.
The “Large” lattice has the same lattice spacing as the “Coarse” lattice but
has a much larger physical size, $L=16.4~{}\mathrm{fm}$. We use the difference
between the “Large” and the “Coarse” results of
$a_{\mu}(R_{\text{max}}<4~{}\mathrm{fm})\times 10^{10}$ as the finite volume
correction. Comparing the results of “Match” and “Coarse”, we estimate the
uncertainty due to the non-zero lattice spacing effects of this model lattice
calculation can be about 12%. Combined with the additional 20% uncertainty due
to the inaccuracy of the model itself (and potential finite volume corrections
of other heavier intermediate states), we obtain our final estimate of the
finite volume correction:

$\displaystyle a_{\mu}^{\text{FV-corr}}(R_{\text{max}}<4~{}\mathrm{fm})\times
10^{10}=-0.47(11)_{\text{syst}}.$ (71)

Similar to the long-distance $\pi^{0}$-exchange contribution, we use the ratio
between the connected and disconnected contributions as described in Eq. 42 to
assign this contribution to the connected and disconnected diagrams. The
results are shown in Tab. 2 as “48I light con FV-corr”, “48I light discon FV-
corr”, and “48I light FV-corr”.

Also, note that we can calculate the long-distance $\pi^{0}$-pole contribution
using this LMD model. In Tab. 4, we also list the contribution to the region
$R_{\text{max}}>4~{}\mathrm{fm}$. Using the LMD model results calculated with
the “Large” lattice, we obtain

$\displaystyle
a_{\mu}^{\pi^{0}\text{-pole};\text{LMD}}(R_{\text{max}}>4~{}\mathrm{fm})\times
10^{10}=2.34.$ (72)

This result agrees with the long-distance $\pi^{0}$-exchange contribution
calculated with $\pi^{0}$ transition form factors from the previously
described lattice QCD calculation and long-distance approximation.

Label | Size | $a^{-1}/\mathrm{GeV}$ | $a_{\mu}(R_{\text{max}}<4~{}\mathrm{fm})\times 10^{10}$ | $a_{\mu}(R_{\text{max}}>4~{}\mathrm{fm})\times 10^{10}$
---|---|---|---|---
Match | $48^{3}\times 96$ | 1.73 | 5.19 | 0.22
Coarse | $24^{3}\times 48$ | 0.865 | 4.65 | 0.20
Large | $64^{3}\times 128$ | 0.865 | 4.18 | 2.34
Table 4:  Lattice calculation of the $\pi^{0}$-pole contribution using the LMD
model. Note that this is not a lattice QCD calculation.

###  IV.5 Pion mass extrapolation

While the calculation is performed very near the physical pion mass, there is
a slight mis-tuning of the lattice parameters which leads to a slightly
heavier $m_{\pi}=139~{}\mathrm{MeV}$ compared to the physical pion mass
$m_{\pi}=134.9766~{}\mathrm{MeV}$. We correct this small difference using the
32D and 24DH ensembles. These two ensembles have the same lattice spacing but
different pion masses. The results of these two ensembles are listed in Tab.
5. Note that for the long-distance contribution
($R_{\text{max}}>4~{}\mathrm{fm}$) for the 32D physical pion mass ensemble, we
use the same long-distance $\pi^{0}$-exchange contribution calculated with
ensemble 48I as described in Section IV.3. Similar to the 48I calculation, we
calculate the “no-pion” contribution and obtain the “hybrid-2fm” results for
the “discon” and “total” results. To reduce the statistical error, we use the
“hybrid-2fm” results to calculate the pion mass correction.

name |  $R_{\text{max}}$ limit | $a_{\mu}(\text{32D})$ | $\times 10^{10}$ | $a_{\mu}(\text{24DH})$ | $\times 10^{10}$
---|---|---|---|---|---
con | 2 fm | $8.17$ | $(0.32)_{\text{stat}}$ | $8.29$ | $(0.10)_{\text{stat}}$
con | 4 fm | $21.29$ | $(3.37)_{\text{stat}}$ | $12.35$ | $(0.24)_{\text{stat}}$
con | $\infty$ | $28.84$ | $(3.39)_{\text{stat}}$ | $12.35$ | $(0.24)_{\text{stat}}$
no-pion | 2 fm | $3.91$ | $(0.63)_{\text{stat}}$ | $4.18$ | $(0.57)_{\text{stat}}$
no-pion | 4 fm | $4.96$ | $(4.46)_{\text{stat}}$ | $6.74$ | $(4.18)_{\text{stat}}$
discon | 2 fm | $-2.10$ | $(0.59)_{\text{stat}}$ | $-1.91$ | $(0.55)_{\text{stat}}$
discon | 4 fm | $-10.70$ | $(3.84)_{\text{stat}}$ | $-2.34$ | $(4.18)_{\text{stat}}$
discon | $\infty$ | $-16.25$ | $(3.84)_{\text{stat}}$ | $-2.34$ | $(4.18)_{\text{stat}}$
discon hybrid-2fm | 4 fm | $-11.74$ | $(2.52)_{\text{stat}}$ | $-4.90$ | $(0.54)_{\text{stat}}$
discon hybrid-2fm | $\infty$ | $-17.30$ | $(2.53)_{\text{stat}}$ | $-4.90$ | $(0.54)_{\text{stat}}$
total | 2 fm | $6.07$ | $(0.66)_{\text{stat}}$ | $6.38$ | $(0.57)_{\text{stat}}$
total | 4 fm | $10.59$ | $(4.97)_{\text{stat}}$ | $10.01$ | $(4.18)_{\text{stat}}$
total | $\infty$ | $12.59$ | $(4.98)_{\text{stat}}$ | $10.01$ | $(4.18)_{\text{stat}}$
total hybrid-2fm | 4 fm | $9.55$ | $(1.12)_{\text{stat}}$ | $7.45$ | $(0.59)_{\text{stat}}$
total hybrid-2fm | $\infty$ | $11.55$ | $(1.12)_{\text{stat}}$ | $7.45$ | $(0.59)_{\text{stat}}$
Table 5: 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 $\pi^{0}$-exchange contribution in the
$R_{\text{max}}>4~{}\mathrm{fm}$ region. For the 24DH ensemble
($m_{\pi}=341~{}\mathrm{MeV}$), we assume the contribution is negligible in
the $R_{\text{max}}>4~{}\mathrm{fm}$ region.

To obtain the pion mass correction to our main 48I results, we assume the
following pion mass dependence for the “no-pion” contribution:

$\displaystyle a_{\mu}^{\text{no-pion}}(m_{\pi})=C_{1}+C_{2}m_{\pi}^{2}.$ (73)

Due to the pion pole piece in the quark-connected diagram, we assume a more
singular pion mass dependence for the “con” contribution:

$\displaystyle a_{\mu}^{\text{con}}(m_{\pi})=C_{1}+\frac{C_{2}}{m_{\pi}^{2}}.$
(74)

The disconnected diagram contribution can be obtained as a combination of the
“con” and “no-pion” contribution as given in Eq. (45)((46)).

The fitting form in Eq. (74)(plus an additional $m_{\pi}^{2}$ term) was used
for the chiral extrapolation of the connected diagram in Ref. [48] and seems
to provide the most accurate description of the connected diagram’s mass
dependence. The correction from the chiral extrapolation is $0.35(7)\times
10^{-10}$ with the above fitting form. In Ref. [26], the $\pi^{0}$-pole
contribution is calculated to behave as $C_{1}+C_{2}\log^{2}(m_{\pi}^{2})$, If
we estimate the pion mass correction with this form instead, the size of the
correction is $0.20(4)\times 10^{-10}$. In Ref. [69], the HLbL is calculated
with Chiral effective theory, we have calculated the pion mass dependence with
Eq. (13) of this reference. The shift of $a_{\mu}$ due to the pion mass change
from the $139~{}\mathrm{MeV}$ used in the lattice calculation to the physical
point $135~{}\mathrm{MeV}$ is $0.25\times 10^{-10}$. Here, we do not include
the correction to the charged pion loop contribution of $a_{\mu}$ due to the
pion mass shift, as the charged pion loop contribution are significantly
reduced compared with the scalar QED estimate due to the pion dressing [70,
71, 72]. We, therefore, estimate a 50% systematic uncertainty to these
corrections to account for possible inaccuracies of the empirical chiral
extrapolations, plus any systematic effects in obtaining the 32D and 24DH
results, including the effects caused by the hybrid method used to calculate
the long-distance part of the disconnected diagrams, the discretization
effects, the finite volume effects, and so on. The final results for the
correction due to the slight mismatch of the pion mass is:

$\displaystyle a_{\mu}^{\text{$m_{\pi}$-corr}}\times 10^{10}$ (75)
$\displaystyle=(a_{\mu}(m_{\pi}=135~{}\mathrm{MeV})-a_{\mu}(m_{\pi}=139~{}\mathrm{MeV}))\times
10^{10}$ $\displaystyle=0.35(0.07)_{\text{stat}}(0.17)_{\text{syst}}.$

We also show the results in Tab. 2 as “48I light con $m_{\pi}$-corr”, “48I
light discon $m_{\pi}$-corr”, and “48I light $m_{\pi}$-corr”.

###  IV.6 Continuum limit extrapolation

In our previous work [24], we used finite-volume lattices for the entire
calculation, including QED. The QEDL scheme was used for the photon
propagators, and sizable discretization effects were observed. We took the
continuum limit with the 48I and 64I ensembles. These two ensembles have
different lattice spacings but all other properties are almost identical. They
have the same gauge and fermion actions, and almost the same pion mass and
volume in physical units. We used the results from these two ensembles to
perform the continuum extrapolation. The continuum limit value is $30(12)\%$
higher compared to the value on the 48I ensemble for the connected piece and
$58(27)\%$ higher for the disconnected piece. The numbers in parentheses are
statistical uncertainties.

In this work, we have mainly used the 48I ensemble to perform the calculation
with the infinite volume, continuum QED weighting function described in Ref.
[46]. Compared with QEDL, the infinite volume QED weighting function reduces
the finite volume effects from power-law to exponential.

One may expect similar discretization effects for these two calculations since
they should be the same in the large volume limit and only differ by finite
volume effects. However, there are two reasons that the discretization errors
in the current work should be much smaller. First, the infinite volume QED
weighting function is calculated in the continuum with a semi-analytic method,
which eliminates the discretization errors in the QEDL weighting function
which is calculated on a discrete lattice. Second, and more importantly, we
apply the subtraction scheme as described in Ref. [46] and repeated in Eq. 12.
In a test calculation of the leptonic light-by-light contribution to the muon
$g-2$, we found about a factor of $5$ reduction in the discretization error
after this subtraction is applied. We expect a similar reduction in
discretization error in the hadronic case, and this is indeed the case in our
calculation of the strange quark connected diagrams.

Figure 9: Strange quark connected part of $a_{\mu}^{\rm HLbL}$ (lower panel).
Continuum limit (upper curve) and finite lattice spacing results on 48I and
64I ensembles. Corresponding summands (upper panel).

The strange quark connected contribution is shown in Fig. 9. Compared with the
light quark, due to the heavier strange quark mass, the statistical noise in
the long-distance region is very small. The overall signal-to-noise ratio is
much better than the light quark case which allows a precise continuum limit
extrapolation. As is shown in the figure, the limit is 7.9(2.6)% higher than
the 48I results. The continuum limit is obtained by performing an
$\mathcal{O}(a^{2})$ extrapolation to $a\to 0$,

$\displaystyle a_{\mu}(a^{2})=C_{1}+C_{2}a^{2}.$ (76)

Compared with the observed $\sim 30$% discretization error in the QEDL
calculation of the light quark connected diagrams, there is indeed a
significant reduction in the relative discretization error with the subtracted
infinite volume weighting function, despite the heavier strange quark mass,
which usually leads to larger discretization errors.

The study of the continuum limit of the strange quark connected diagrams
strongly suggests the light quark connected diagrams, in particular in the
relatively short distance region where we performed a reliable continuum
extrapolation for the strange quark connected diagrams, is under 8%. In the
long distance (large $R_{\text{max}}$) region, the contributions are mostly
coming from the $\pi^{0}$-exchange contribution. As can be seen from Eq.
(IV.3), the hadronic four-point function in the HLbL diagrams in the long-
distance part can be approximated by the product of two such transition form
factors. Therefore, the discretization effects of the direct calculation of
the four-point function in the long distance region will be similar to the
discretization effects of the calculation of the square of the transition form
factor. In the $\pi^{0}\to e^{+}e^{-}$ calculation of Ref. [73], this relevant
$\pi^{0}\to\gamma\gamma$ transition form factor is studied with the same 48I
ensemble as this calculation. In addition, the continuum limit is calculated
with a parallel 64I ensemble calculation. We quote the corrections obtained in
the calculation in Tab. 6. Note that the correction is zero consistent given
the statistical uncertainty. In Fig. 4 of Ref. [73], the partial sum of these
quantities for both 48I and 64I are plotted. More statistically precise
agreement is observed at smaller time separation of the two electromagnetic
vector currents.

Observable | Percentage correction
---|---
$(\mathrm{Im}\mathcal{A})^{2}$ | -11.4(11.6)%
$(\mathrm{Re}\mathcal{A})^{2}$ | -16.8(15.8)%
Table 6:  Percent correction to obtain the continuum limit relative to the
value computed on the 48I ensemble in the $\pi^{0}\to e^{+}e^{-}$ calculation
in Ref. [73]. The values in parentheses are the statistical uncertainty of the
corresponding quantity. The square of the amplitude is used to estimate the
lattice spacing error of the long-distance $\pi^{0}$-exchange contribution to
the muon $g-2$ in this work.

For the light quark disconnected diagrams, we again separately discuss the
short and long distance regions. In Fig. 5, we observe that in the short
distance region ($R_{\text{max}}\lesssim 2~{}\mathrm{fm}$), the disconnected
diagram contribution is much smaller than the connected contribution. This is
expected as the disconnected diagrams involve two quark loops and the exchange
of at least two gluons. Therefore, the disconnected diagrams are suppressed by
$\mathcal{O}(\alpha_{s}^{2})$ and $1/N_{c}$ relative to the connected
diagrams. Due to the smallness of the contributions from the short-distance
disconnected diagrams, we expect that the same reasoning applies to the
discretization effects.

For the long distance region, we again have Eq. (42), which mandates the
relationship with the connected diagrams, even for non-zero lattice spacing
and finite volume (so long that the lattice is large enough to contain the
correlation function). Therefore, we expect the long distance part of the
disconnected diagrams to have the same relative discretization effects as the
connected diagrams and also the combined total. Note that the situation is
different from the Mainz group’s work in Ref. [47, 48]. While the argument
presented here should also apply to the hadronic function in their
calculation, the different QED kernels (due to the subtraction scheme of the
QED kernel and the choice of $x_{\text{ref}}$) used for the connected and
disconnected diagrams lead to different shapes of the summands and also
different discretization errors.

Without a 64I ensemble calculation of similar precision, the continuum limit
cannot be reliably taken. However, based on the evidence presented above, in
particular 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 strange quark disconnected diagrams. The corrections
for light quark are shown in Tab. 2 as “48I light con $a^{2}$-corr”, “48I
light discon $a^{2}$-corr”, and “48I light $a^{2}$-corr”. The corrections for
the strange quark are shown in Tab. 3 as “48I strange con $a^{2}$-corr”, “48I
strange discon $a^{2}$-corr”, and “48I strange $a^{2}$-corr”.

###  IV.7 Sub-leading disconnected diagrams

So far, our discussions have centered on the first two diagrams in Fig. 2. The
remaining diagrams, which are expected to be small due to the suppression by
flavor SU(3) and quark electric charge factors, are not yet included. In our
earlier work [24], we have already calculated the third diagram in Fig. 2,
which is expected to be the largest within the remaining the sub-leading
disconnected diagrams based on the flavor SU(3) and quark electric charge
factor counting. To estimate the contribution from this diagram, we employ the
same subtracted infinite volume QED weighting function and the 24D ensemble
for the quark loops. The result, shown in Fig. 10, is zero consistent, and we
estimate this contribution to be $0.0(0.5)\times 10^{-10}$.

Figure 10: Sub-leading disconnected diagrams computed on the 24D ensemble.
[24]

In the more recent work by the Mainz group [48], all the sub-leading diagrams
have been calculated. The result is still zero consistent but a more stringent
bound is obtained. Therefore, in this work, we use their result to account for
the contribution from all the sub-leading disconnected diagrams. The value is
shown in Tab. 7 as “sub-leading discon”.

###  IV.8 Charm quark contributions

The charm quark contribution is also expected to be very small due its heavy
mass relative to the light quark and the strange quark. In the 2020 white
paper [5], its contribution is included as $0.3(0.1)\times 10^{-10}$,
calculated 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 mass of the charm quark and larger electric
charge. In a recent work by the Mainz group [49], the charm quark contribution
is calculated with lattice QCD, including both the charm quark connected and
disconnected diagrams. The result for the total charm quark contribution is
$0.28(0.05)\times 10^{-10}$, where the uncertainty is mostly due to the
systematic effects from modeling the lattice spacing and $M_{\eta_{c}}$
dependence. In this work, we use this more recent lattice calculated result to
account for the contribution from the charm quark. The values are shown in
Tab. 7 as “charm con”, “charm discon” “charm”.

Contribution | $a_{\mu}$ | $\times 10^{10}$
---|---|---
light con | $25.70$ | $(1.33)_{\text{stat}}(1.99)_{\text{syst}}~{}[2.39]$
light discon | $-12.71$ | $(1.87)_{\text{stat}}(1.17)_{\text{syst}}~{}[2.20]$
light discon hybrid-2.5fm | $-13.50$ | $(1.36)_{\text{stat}}(1.21)_{\text{syst}}~{}[1.82]$
light discon hybrid-2fm | $-13.37$ | $(1.29)_{\text{stat}}(1.46)_{\text{syst}}~{}[1.95]$
light total | $12.99$ | $(2.11)_{\text{stat}}(0.90)_{\text{syst}}~{}[2.29]$
light total hybrid-2.5fm | $12.20$ | $(1.01)_{\text{stat}}(0.95)_{\text{syst}}~{}[1.38]$
light total hybrid-2fm | $12.33$ | $(0.90)_{\text{stat}}(1.25)_{\text{syst}}~{}[1.55]$
strange con | $0.35$ | $(0.01)_{\text{stat}}$
strange discon | $-0.38$ | $(0.61)_{\text{stat}}(0.03)_{\text{syst}}~{}[0.61]$
strange discon hybrid-2.5fm | $-0.36$ | $(0.22)_{\text{stat}}(0.03)_{\text{syst}}~{}[0.22]$
strange discon hybrid-2fm | $-0.28$ | $(0.12)_{\text{stat}}(0.04)_{\text{syst}}~{}[0.13]$
strange total | $-0.03$ | $(0.61)_{\text{stat}}(0.03)_{\text{syst}}~{}[0.61]$
strange total hybrid-2.5fm | $-0.00$ | $(0.22)_{\text{stat}}(0.03)_{\text{syst}}~{}[0.23]$
strange total hybrid-2fm | $0.07$ | $(0.12)_{\text{stat}}(0.04)_{\text{syst}}~{}[0.13]$
sub-leading discon | $0.00$ |  $\phantom{(0.00)_{\text{stat}}}(0.07)_{\text{syst}}$ [48]
charm con | $0.31$ |  $\phantom{(0.00)_{\text{stat}}}(0.04)_{\text{syst}}$ [49]
charm discon | $-0.03$ |  $\phantom{(0.00)_{\text{stat}}}(0.02)_{\text{syst}}$ [49]
charm total | $0.28$ |  $\phantom{(0.00)_{\text{stat}}}(0.05)_{\text{syst}}$ [49]
con | $26.36$ | $(1.33)_{\text{stat}}(1.99)_{\text{syst}}~{}[2.39]$
discon | $-13.12$ | $(2.30)_{\text{stat}}(1.18)_{\text{syst}}~{}[2.59]$
discon hybrid-2.5fm | $-13.89$ | $(1.47)_{\text{stat}}(1.22)_{\text{syst}}~{}[1.91]$
discon hybrid-2fm | $-13.68$ | $(1.35)_{\text{stat}}(1.47)_{\text{syst}}~{}[1.99]$
total | $13.24$ | $(2.53)_{\text{stat}}(0.90)_{\text{syst}}~{}[2.68]$
total hybrid-2.5fm | $12.47$ | $(1.15)_{\text{stat}}(0.95)_{\text{syst}}~{}[1.49]$
total hybrid-2fm | $12.68$ | $(0.98)_{\text{stat}}(1.26)_{\text{syst}}~{}[1.59]$
Table 7: Summary of final results. Values are different contributions to
$a_{\mu}\times 10^{10}$ where $a_{\mu}=(g_{\mu}-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.1 and IV.2. The result for “total hybrid-2.5fm” is
used as our final result.

##  V Conclusion

We summarized all the results discussed above in Tabs. 2, 3, and 7. Adding all
the individual contributions and corrections, we obtain our final result for
the HLbL contribution to the muon $g-2$:

$\displaystyle a_{\mu}^{\text{HLbL}}\times
10^{10}=12.47(1.15)_{\text{stat}}(0.95)_{\text{syst}}~{}[1.49],$ (77)

where the systematic errors were added in quadrature. We also give the HLbL
contributions to the muon $g-2$ from the connected and disconnected pieces
separately:

$\displaystyle a_{\mu}^{\text{HLbL,con}}\times 10^{10}$
$\displaystyle=26.36(1.33)_{\text{stat}}(1.99)_{\text{syst}}~{}[2.39],$ (78)
$\displaystyle a_{\mu}^{\text{HLbL,discon}}\times 10^{10}$
$\displaystyle=-13.89(1.47)_{\text{stat}}(1.22)_{\text{syst}}~{}[1.91].$ (79)

Numbers in square square brackets denote total erorrs, combining the
statistical uncertainty (“stat”) and systematic ones (“sys”) in quadrature.

Here, we summarize techniques used in this calculation that we believe are
important to obtain the precision for HLbL scattering at the physical pion
mass:

  * •

The subtracted infinite volume QED weighting function developed in our
previous work [46].

  * •

Use all combinations of the point source propagators to calculate the HLbL
diagrams as described in section III.

  * •

The adaptive sampling scheme used for the disconnected diagram calculation as
described in section III.

  * •

The rearrangement of the connected and disconnected diagrams in Eq. (43,46)
based on Eq. (42).

  * •

The AMA method [58] and efficient GPU solver for MDWF propagators from Grid
and GPT.

| $a_{\mu}$ | $\times 10^{10}$
---|---|---
all con QEDL | $24.46$ | $(2.35)_{\text{stat}}(5.11)_{\text{syst}}~{}[5.62]$
all con diff | $1.90$ | $(2.76)_{\text{stat}}(5.48)_{\text{syst}}~{}[6.14]$
all discon QEDL | $-16.45$ | $(2.09)_{\text{stat}}(3.99)_{\text{syst}}~{}[4.50]$
all discon diff | $2.56$ | $(2.57)_{\text{stat}}(4.17)_{\text{syst}}~{}[4.90]$
total QEDL | $8.17$ | $(3.03)_{\text{stat}}(1.77)_{\text{syst}}~{}[3.51]$
total diff | $4.30$ | $(3.25)_{\text{stat}}(2.01)_{\text{syst}}~{}[3.82]$
Table 8: Comparison of this work with our previous QEDL results [24] plus the
charm quark contribution from [5, 23]. The field “all con” includes all
connected diagram contributions, including the charm quark. The field “all
discon” includes all the disconnected diagram contributions, including the
sub-leading disconnected diagrams and disconnected diagrams including charm
quarks. The fields “all con diff”, “all discon diff”, “total diff” shows the
difference obtained by subtracting the QEDL results from the new results
presented in this work. The statistical and systematic uncertainties of the
two works are almost independent. We therefore add them in quadrature for the
difference. The numbers in the square bracket are the statistical and
systematic uncertainty combined in quadrature.  Figure 11: Comparison of our
result with values in the literature. The hadronic model values are from Refs.
[26, 27, 28, 29]. The dispersive data driven result is compiled in Ref. [5].
The lattice 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. The comparison is summarized in Table 8. We can see
that the results for both the connected and disconnected diagrams are in good
agreement. For the total, the current result is $1.12$ standard deviations
higher than the previous results, possibly due to a slightly larger
statistical fluctuation. We also compare the final result in this work with
the existing literature in Fig. 11. The new result is consistent with previous
determinations.

###### Acknowledgements.

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 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-SC0011941.
M.H. is supported by Japan Grants-in-Aid for Scientific Research, No20K03926.
T.I., C.J., and C.L. were supported in part by US DOE Contract
DESC0012704(BNL). T.I. and C.J. were supported in part by the Scientific
Discovery through Advanced Computing (SciDAC) program LAB 22-2580. T.I. is
also supported by the Department of Energy, Laboratory Directed Research and
Development (LDRD No. 23-051) of BNL and RIKEN BNL Research Center, and by
JSPS KAKENHI under grant numbers JP26400261, JP17H02906. C.L. has been
supported by a DOE Office of Science Early Career Award. We developed the
computational code used for this work based on the Columbia Physics System
(CPS), Grid, GPT, and Qlattice. Computations were performed mainly under the
ALCC Program of the US DOE on the SUMMIT computer at the Oak Ridge Leadership
Computing Facility, a DOE Office of Science Facility supported under Contract
No. DE-AC05-00OR22725. These calculations used gauge configurations and
propagators created using resources of the Argonne Leadership Computing
Facility, which is a DOE Office of Science User Facility supported under
Contract DE-AC02-06CH11357. We also acknowledge computer resources at the
Oakforest-PACS supercomputer system at Tokyo University, partly through the
HPCI System Research Project (hp180151, hp190137), the BNL SDCC computer
clusters at Brookhaven National Lab as well as computing resources provided
through USQCD at Brookhaven and Jefferson National Labs.

##  Appendix A Notation

We use $S_{{\mu}}$ and $G$ to denote free muon and photon propagators:

$\displaystyle S_{{\mu}}(x,y)$
$\displaystyle=\int\frac{d^{4}p}{(2\pi)^{4}}\frac{1}{i\\!\not{\\!p}+m}e^{ip\cdot(x-y)}$
(80)
$\displaystyle=\left(-\\!\not{\\!\partial}_{x}+m\right)\int\frac{d^{4}p}{(2\pi)^{4}}\frac{1}{p^{2}+m^{2}}e^{ip\cdot(x-y)},$
(81) $\displaystyle G(x,y)$
$\displaystyle=\int\frac{d^{4}p}{(2\pi)^{4}}\frac{1}{p^{2}}e^{ip\cdot(x-y)}$
(82) $\displaystyle=\frac{1}{4\pi^{2}}\frac{1}{(x-y)^{2}}.$ (83)

The $\gamma_{\mu}$ matrices satisfy the Euclidean space-time metric

$\displaystyle\gamma_{\mu}\gamma_{\nu}+\gamma_{\nu}\gamma_{\mu}=2\delta_{\mu,\nu}\,.$
(84)

##  Appendix B $\pi^{0}$ long-distance contribution

In Eq. (IV.3), we replaced the QCD, Euclidean space-time, four-current
connected Green’s function with the product of two amplitudes, each coupling a
pair of currents to an on-shell $\pi^{0}$. These two amplitudes are joined by
a pion propagator and all amplitudes are expressed in position space, so they
can be directly inserted in our standard position-space evaluation of the HLbL
amplitude. Since the final expression involves two independent factors
evaluating the $\pi\gamma\gamma$ coupling which are connected by an analytic,
position-space pion propagator, this QCD part of the HLbL amplitude can be
evaluated in a “QCD volume” of arbitrary size. In particular, this volume
could be much larger than that of the gauge configurations used to compute
each $\pi\gamma\gamma$ vertex function. Here we work out a concrete derivation
of this formula that can be used to evaluate the long-distance part of the
$\pi^{0}$ exchange contribution to leading order in $1/L$. We leave open the
possibility that this approach could be developed further to systematically
capture terms falling with higher powers of $1/L$ if the large volume
$\pi^{0}$ contribution is expressed as a power series in $1/L^{n}$ where $L$
is the size of the QCD volume. We assume the QED volume to be infinite.

We begin with the Euclidean-space Green’s function defined in Eq. (II):

$\displaystyle
6e^{4}\mathcal{H}_{\mu^{\prime},\mu,\nu^{\prime},\nu}(x^{\prime},x,y^{\prime},y)$
(85)
$\displaystyle=\bigl{\langle}T\bigl{(}J_{\mu^{\prime}}(x^{\prime})J_{\mu}(x)J_{\nu^{\prime}}(y^{\prime})J_{\nu}(y)\bigr{)}\bigr{\rangle}.$

We will choose $x^{\prime}$ and $x$ close to each other as are $y^{\prime}$
and $y$. However, we will assume that the $(x^{\prime},x)$ pair is far from
the $(y^{\prime},y)$ pair. Define:

$\displaystyle\widetilde{x}$ $\displaystyle=x^{\prime}-x,$ (86)
$\displaystyle\widetilde{y}$ $\displaystyle=y^{\prime}-y.$ (87)

Since Euclidean space-time is rotational invariant, without loss of
generality, we can choose $x-y$ to be along the Euclidean time direction, with
$x-y=(t,\vec{0})$ and $t>0$.

We now follow the usual steps to obtain a variant of the Källen-Lehman
representation of the amplitude but keep only the single $\pi^{0}$
intermediate state since, as the lightest particle, its exchange will dominate
this Green’s function when $x$ and $y$ are far separated:

$\displaystyle
6e^{4}\mathcal{H}^{\pi^{0}}_{\mu^{\prime},\mu,\nu^{\prime},\nu}(x^{\prime},x,y^{\prime},y)$
(88) $\displaystyle=$
$\displaystyle\int\\!\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{2E_{\pi,\vec{p}}}\bigl{\langle}0\bigl{|}T\bigl{(}J_{\mu^{\prime}}(x^{\prime})J_{\mu}(x)\bigr{)}\bigr{|}\pi^{0}(\vec{p})\bigr{\rangle}$
$\displaystyle\hskip
56.9055pt\times\bigl{\langle}\pi^{0}(\vec{p})\bigl{|}T\bigl{(}J_{\nu^{\prime}}(y^{\prime})J_{\nu}(y)\bigr{)}\bigr{|}0\bigr{\rangle},$

where the superscript $\pi^{0}$ indicates that only the $\pi^{0}$ contribution
to $\mathcal{H}$ is represented.

Next we use four-dimensional translational invariance to remove the variables
$x$ and $y$ from the two current-current-$\pi^{0}$ amplitudes. We will also
replace these two $O(4)$-covariant current-current-$\pi^{0}$ amplitudes by
functions of the four momentum $p=(iE_{\pi,\vec{p}},\vec{p})$. Given the
factors of $\exp{(ip\cdot x)}$ and $\exp{(-ip\cdot y)}$ which appear, we can
then replace the four vector $p_{\mu}$ by $-i\partial/\partial{x}_{\mu}$ or
$i\partial/\partial{y}_{\mu}$ as needed. For example,

$\displaystyle\bigl{\langle}0\bigl{|}T\bigl{(}J_{\mu^{\prime}}(x^{\prime})J_{\mu}(x)\bigr{)}\bigr{|}\pi^{0}(\vec{p})\bigr{\rangle}$
(89)
$\displaystyle\hskip-108.405pt=\bigl{\langle}0\bigl{|}T\bigl{(}J_{\mu^{\prime}}(\widetilde{x})J_{\mu}(0)\bigr{)}\bigr{|}\pi^{0}(\vec{p})\bigr{\rangle}e^{ip\cdot
x}$
$\displaystyle\hskip-108.405pt=\mathcal{F}_{\mu^{\prime},\mu}(\widetilde{x},p)e^{ip\cdot
x}$
$\displaystyle\hskip-108.405pt=\mathcal{F}_{\mu^{\prime},\mu}\left(\widetilde{x},-i\frac{\partial}{\partial{x}}\right)e^{ip\cdot
x}.$

where the definition of the $\pi^{0}$ transition form factors follows Eq.
(47).

Using this approach to remove the explicit dependence of the two amplitudes on
$p_{\mu}$, we can then perform the final step of the usual Källen-Lehman
derivation and introduce the free pion propagator:

$\displaystyle\int\\!\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{2E_{\pi,\vec{p}}}e^{ip\cdot(x-y)}$
$\displaystyle=\int\\!\frac{d^{3}p}{(2\pi)^{3}}\frac{1}{2E_{\pi,\vec{p}}}e^{i\vec{p}\cdot(\vec{x}-\vec{y})}e^{-E_{\pi,\vec{p}}({x}_{t}-{y}_{t})}$
(90)
$\displaystyle=\int\\!\frac{d^{4}p}{(2\pi)^{4}}\frac{e^{ip\cdot(x-y)}}{p^{2}+m_{\pi}^{2}}$
(91) $\displaystyle=D_{\pi^{0}}(x-y)$ (92)

where $D_{\pi^{0}}(x-y)$ is the free Euclidean-space propagator for a scalar
particle of mass $m_{\pi}$.

We can then rewrite Eq. (88) in terms of $\mathcal{F}$ and $D_{\pi^{0}}$ to
obtain the result:

$\displaystyle
6e^{4}\mathcal{H}^{\pi^{0}}_{\mu^{\prime},\mu,\nu^{\prime},\nu}(x^{\prime},x,y^{\prime},y)$
(93)
$\displaystyle=\mathcal{F}_{\mu^{\prime},\mu}\left(\widetilde{x},-i\frac{\partial}{\partial{x}_{\mu}}\right)\mathcal{F}_{\nu^{\prime},\nu}\left(\widetilde{y},i\frac{\partial}{\partial{y}_{\mu}}\right)D_{\pi^{0}}(x-y).$

At this step, we have made explicit the $O(4)$-covariance of the right-hand
side of the equation. So this equation should hold for an arbitrary
orientation of $x-y$.

So far all of the steps taken have been exact. We expect the quantity
computed, the contribution of a single pion exchange, to dominate the long
distance limit in which $|x-y|$ is large, i.e. $|x-y|\gtrsim L$. Now we will
evaluate the derivatives with respect to $x$ and $y$ which appear in these
equations but keep only the leading term in an expansion in powers of $1/L$.
For example, to leading order in $1/|x-y|$:

$\displaystyle\quad\prod_{i=1}^{N}\left(\frac{\partial}{\partial{x}_{\rho_{i}}}\right)D_{\pi^{0}}(x-y)$
$\displaystyle=\prod_{i=1}^{N}\left(\frac{\partial}{\partial{x}_{\rho_{i}}}\right)\frac{c\,e^{-m_{\pi}|x-y|}}{|x-y|^{3/2}}$
(94)
$\displaystyle\approx\prod_{i=1}^{N}\left(-m_{\pi}\frac{(x-y)_{\rho_{i}}}{|x-y|}\right)\frac{c\,e^{-m_{\pi}|x-y|}}{|x-y|^{3/2}},$
(95)

where the relation holds for arbitrary $N$ and $\rho_{i}$ labels the component
of $x$ appearing in the $i^{th}$ derivative in the product of $N$ derivatives.
Here we use the large distance behavior of $D_{\pi^{0}}(x-y)$:

$\displaystyle
D_{\pi^{0}}(x-y)\approx\frac{c\,e^{-m_{\pi}|x-y|}}{|x-y|^{3/2}}.$ (96)

Introducing this approximation into Eqs. (93), we obtain

$\displaystyle
6e^{4}\mathcal{H}^{\pi^{0}}_{\mu^{\prime},\mu,\nu^{\prime},\nu}(x^{\prime},x,y^{\prime},y)$
(97)
$\displaystyle\approx\mathcal{F}_{\mu^{\prime},\mu}\left(\widetilde{x},im_{\pi}\hat{n}\right)\mathcal{F}_{\nu^{\prime},\nu}\left(\widetilde{y},-im_{\pi}\hat{n}\right)D_{\pi^{0}}(x-y).$

where $\hat{n}_{\mu}$ is the Euclidean unit vector $(x-y)_{\mu}/|x-y|$. This
is finally the equation of interest, Eq. (IV.3), which appeared earlier. Since
the two quantities expressed in terms of the form factor $\mathcal{F}$ can be
evaluated from independent lattice QCD ensemble averages, the expression on
the right hand side of Eq. (97) can be evaluated from a standard lattice QCD
ensemble average followed by an $O(4)$ rotation as defined in Eq. (57), which
we repeat here:

$\displaystyle\mathcal{F}_{\mu^{\prime},\mu}\Big{(}\widetilde{x},im_{\pi}\hat{n}\Big{)}=\Lambda_{\mu^{\prime},\rho^{\prime}}\Lambda_{\mu,\rho}\mathcal{F}_{\rho^{\prime},\rho}\Big{(}\widetilde{x}^{\prime},im_{\pi}\hat{t}\Big{)},$
(98)

The $O(4)$ rotation matrix $\Lambda=\Lambda(\hat{n})$ and Euclidean space-time
coordinate $\tilde{x}^{\prime}$ satisfy:

$\displaystyle\delta_{\mu,\nu}$
$\displaystyle=\Lambda_{\mu,\mu^{\prime}}\Lambda_{\nu,\nu^{\prime}}\delta_{\mu^{\prime},\nu^{\prime}},$
(99) $\displaystyle\hat{n}_{\mu}$
$\displaystyle=\Lambda_{\mu,\mu^{\prime}}\hat{t}_{\mu^{\prime}},$ (100)
$\displaystyle\widetilde{x}_{\mu}$
$\displaystyle=\Lambda_{\mu,\mu^{\prime}}\widetilde{x}^{\prime}_{\mu^{\prime}}.$
(101)

The rotated form factor
$\mathcal{F}_{\rho^{\prime},\rho}\Big{(}\widetilde{x}^{\prime},im_{\pi}\hat{t}\Big{)}$
follows the definition Eq. (89) or Eq. (47) with zero momentum $\pi^{0}$
state:

$\displaystyle\vspace{-1.0cm}\mathcal{F}_{\rho^{\prime},\rho}\Big{(}\widetilde{x}^{\prime},im_{\pi}\hat{t}\Big{)}$
$\displaystyle=\langle
0|TJ_{\rho^{\prime}}(\widetilde{x}^{\prime})J_{\rho}(0)|\pi^{0}(\vec{p}=\vec{0})\rangle$
$\displaystyle=\sqrt{2m_{\pi}L^{3}}\lim_{t_{\text{tsep}}\to+\infty}\frac{\langle
TJ_{\rho^{\prime}}(\widetilde{x}^{\prime})J_{\rho}(0)\pi^{0}(-t_{\text{sep}})\rangle}{\sqrt{\langle
T\pi^{0}(t_{\text{sep}})\pi^{0}(-t_{\text{sep}})\rangle}},$ (102)

where $\pi^{0}(t)$ is a zero momentum $\pi^{0}$ interpolating operator. We use
Coulomb fixed gauge wall source pion operator in our numerical lattice
calculation. The other factor
$\mathcal{F}_{\nu^{\prime},\nu}\left(\widetilde{y},-im_{\pi}\hat{n}\right)$
can be calculated similarly.

The above procedure allows a direct calculation of the $\pi^{0}$ contribution
to the long distance part the HLbL amplitude without relying on any
parameterization or modelling of the $\pi^{0}$ transition form factor.

In this calculation, we employ this long-distance $\pi^{0}$-exchange
approximation in the following region:

$\displaystyle
R_{\text{max}}=\text{max}(|x-y|,|y^{\prime}-y|,|x-y^{\prime}|)\geq
4~{}\mathrm{fm},$ (103)

as defined in Eq. (41). This long-distance hadronic function is then combined
with infinite volume QED weighting function for HLbL as illustrated in Fig. 8.
Note that we need to sum over all permutations of $x,y,y^{\prime}$.

##  Appendix C Summand fitting form

In section IV.1, we have combined the disconnected diagrams and connected
diagrams with a certain fraction multiplied to form $a_{\mu}^{\text{no-pion}}$
according to Eq. (43). We have fitted the long distance behavior of this
combination with Eq. (44) and the results are shown in Fig. 6. In this
appendix, we will explore other possible fitting forms and will also apply
these fitting functions to other contributions.

The fitting forms are

  * •

fit1: This is the form used in the main calculation, same as Eq. (43).

$\displaystyle
f(R_{\text{max}})=A\frac{R_{\text{max}}^{6}}{R_{\text{max}}^{3}+C^{3}}e^{-BR_{\text{max}}}$
(104)

  * •

fit2:

$\displaystyle f(R_{\text{max}})=AR_{\text{max}}^{3}e^{-BR_{\text{max}}}$
(105)

  * •

fit3:

$\displaystyle
f(R_{\text{max}})=A\frac{R_{\text{max}}^{6}}{(R_{\text{max}}^{3}+C^{3})^{2}}e^{-BR_{\text{max}}}$
(106)

  * •

fit4:

$\displaystyle
f(R_{\text{max}})=A\frac{R_{\text{max}}^{6}}{R_{\text{max}}^{3}+C^{3}}e^{-m_{\eta}R_{\text{max}}}$
(107)

  * •

fit5:

$\displaystyle
f(R_{\text{max}})=AR_{\text{max}}^{3}e^{-m_{\eta}R_{\text{max}}}$ (108)

  * •

fit6:

$\displaystyle
f(R_{\text{max}})=A\frac{R_{\text{max}}^{6}}{(R_{\text{max}}^{3}+C^{3})^{2}}e^{-m_{\eta}R_{\text{max}}}$
(109)

  * •

fit7:

$\displaystyle
f(R_{\text{max}})=A\frac{R_{\text{max}}^{6}}{(R_{\text{max}}^{3}+C^{3})^{2}}e^{-2m_{\pi}R_{\text{max}}}$
(110)

  * •

fit8:

$\displaystyle
f(R_{\text{max}})=A\frac{R_{\text{max}}^{6}}{(R_{\text{max}}^{3}+C^{3})^{2}}e^{-m_{\pi}R_{\text{max}}}$
(111)

where we set $m_{\pi}=139~{}\mathrm{MeV}$ and $m_{\eta}=550~{}\mathrm{MeV}$.
For “fit1”, “fit2”, “fit3”, we always constrain the $B$ parameter to be larger
or equal to $m_{\pi}$.

In Figs. 12,13,14,15,16,17, and 18, we show the resulting fits along with the
data.

We also apply these fit functions to different contributions, including the
strange quark connected diagram, the light quark connected diagrams only, and
the light quark disconnected diagrams only. The strange quark connected fits
are displayed in Figs. 19,20,21 and, for the light quark connected diagrams,
in Figs. 22,23,24,25. And finally, the light quark disconnected results are
shown in Figs. 26,27,28,29. The first fit function (fit1) defined in Eq. (104)
fits the contribution from the strange quark connected diagrams, which are
statistically very precise, almost perfectly. This is the major reason that
this fit function is used as central value for the tail part of the light
quark no-pion contribution in our main calculation described in section IV.1.

Table 9 summarizes the fits. For the light quark no-pion contributions, we
observe roughly 100% variation about the central value (“fit1”), which is the
reason for the 100% systematic uncertainty assigned to the contributions “48I
light no-pion $R_{\text{max}}>2.5\mathrm{fm}$” and “48I light no-pion
$R_{\text{max}}>2\mathrm{fm}$” in Table 2. Note that the “fit3” result for
“48I light no-pion” has large statistical uncertainty. This is due to the fact
that the data does not constrain the exponent $B$ very much, in contrast to
“fit6” and “fit7” which describe the data well. Since the intermediate states
should have at least energy $2m_{\pi}$, we view the results from “fit7” to be
the upper bound for the fitting form “fit3”, which falls into the range of the
results obtained with “fit1” and the 100% systematic uncertainty based on the
other fits.

From Tab. 9, we observe that the fits to the long distance parts of “48I light
con” and “48I light discon” are very unreliable. Note that the statistical
error for the connected diagrams and disconnected diagrams are almost
independent. If we fit the contributions from the connected and disconnected
diagrams individually and add the results, the statistical error will
approximately add up in quadrature. The difficulty of fitting the
contributions from the individual connected and disconnected diagrams is
likely due to the long tail from $\pi^{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. (43) cancels the long distance $\pi^{0}$ exchange
contribution entirely, making the fit much more reliable. Still, we found some
fitting form dependence even after the cancellation of the $\pi^{0}$ exchange
contribution, which we accommodate with the 100% systematic uncertainty
estimated previously.

Contribution name | Fitting form | $R_{\text{max}}>4\mathrm{fm}$ | $R_{\text{max}}>2.5\mathrm{fm}$ | $R_{\text{max}}>2.0\mathrm{fm}$
---|---|---|---|---
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 $\times 1000$ | fit1 | $0.06(0.00)$ | $5.90(0.06)$ | $22.72(0.15)$
48I strange con $\times 1000$ | fit2 | $0.08(0.00)$ | $6.13(0.05)$ | $22.64(0.14)$
48I strange con $\times 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.45)$ | $-8.73(3.61)$ | $-10.78(3.67)$
48I light discon | fit8 | $-5.22(1.84)$ | $-11.46(3.08)$ | $-13.46(3.24)$
Table 9: The long distance part of the $a_{\mu}^{\text{no-pion}}\times
10^{10}$ as defined in Eq. 43 from different fits to the 48I ensemble data.
The numbers in parentheses are statistical uncertainties. 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.

Figure 12: Same as Fig. 6, reproduced for easy comparison with the other
figures in this appendix. The fit function is Eq. (104). The fit starts at
$0.5~{}\mathrm{fm}$.

Figure 13: Similar to Fig. 6 but with the fit function in Eq. (105). The fit
starts at $1.0~{}\mathrm{fm}$.

Figure 14: Similar to Fig. 6 but with the fit function in Eq. (106). The fit
starts at $0.5~{}\mathrm{fm}$.

Figure 15: Similar to Fig. 6 but with the fit function in Eq. (107). The fit
starts at $0.5~{}\mathrm{fm}$.

Figure 16: Similar to Fig. 6 but with the fit function in Eq. (108). The fit
starts at $1.0~{}\mathrm{fm}$.

Figure 17: Similar to Fig. 6 but with the fit function in Eq. (109). The fit
starts at $0.5~{}\mathrm{fm}$.

Figure 18: Similar to Fig. 6 but with the fit function in Eq. (110). The fit
starts at $0.5~{}\mathrm{fm}$.

Figure 19: Similar to Fig. 6 but for the strange quark connected contriubtion
with the fit function in Eq. (104). The fit starts at $0.5~{}\mathrm{fm}$.

Figure 20: Similar to Fig. 6 but for the strange quark connected contribution
with the fit function in Eq. (105). The fit starts at $1.0~{}\mathrm{fm}$.

Figure 21: Similar to Fig. 6 but for the strange quark connected contribution
with the fit function in Eq. (106). The fit starts at $0.5~{}\mathrm{fm}$.

Figure 22: Similar to Fig. 6 but for light quark connected contribution with
the fit function in Eq. (104). The fit starts at $1.5~{}\mathrm{fm}$.

Figure 23: Similar to Fig. 6 but for light quark connected contribution with a
the fit function in Eq. (105). The fit starts at $1.5~{}\mathrm{fm}$.

Figure 24: Similar to Fig. 6 but for the light quark connected contribution
with the fit function in Eq. (106). The fit starts at $1.5~{}\mathrm{fm}$.

Figure 25: Similar to Fig. 6 but for the light quark connected contribution
with the fit function in Eq. (111). The fit starts at $1.5~{}\mathrm{fm}$.

Figure 26: Similar to Fig. 6 but for the light quark disconnected contribution
with the fit function in Eq. (104). The fit starts at $1.5~{}\mathrm{fm}$.

Figure 27: Similar to Fig. 6 but for the light quark disconnected conrtibution
with the fit function in Eq. (105). The fit starts at $1.5~{}\mathrm{fm}$.

Figure 28: Similar to Fig. 6 but for the light quark disconnected contribution
with the fit function in Eq. (106). The fit starts at $1.5~{}\mathrm{fm}$.

Figure 29: Similar to Fig. 6 but for the light quark disconnected contribution
with the fit function in Eq. (111). The fit starts at $1.5~{}\mathrm{fm}$.

## References

  * Abi  _et al._ [2021] B. Abi _et al._ (Muon g-2), Phys. Rev. Lett. 126, 141801 (2021), arXiv:2104.03281 [hep-ex] .
  * Abe  _et al._ [2019] M. Abe _et al._ , PTEP 2019, 053C02 (2019), arXiv:1901.03047 [physics.ins-det] .
  * Sato [2021] Y. Sato (J-PARC E34), JPS Conf. Proc. 33, 011110 (2021).
  * Bennett  _et al._ [2006] G. Bennett _et al._ (Muon G-2), Phys.Rev. D73, 072003 (2006), arXiv:hep-ex/0602035 [hep-ex] .
  * Aoyama  _et al._ [2020] T. Aoyama _et al._ , Phys. Rept. 887, 1 (2020), arXiv:2006.04822 [hep-ph] .
  * Aoyama  _et al._ [2012] T. Aoyama, M. Hayakawa, T. Kinoshita, and M. Nio, Phys. Rev. Lett. 109, 111808 (2012), arXiv:1205.5370 [hep-ph] .
  * Aoyama  _et al._ [2019] T. Aoyama, T. Kinoshita, and M. Nio, Atoms 7, 28 (2019).
  * Czarnecki  _et al._ [2003] A. Czarnecki, W. J. Marciano, and A. Vainshtein, Phys. Rev. D67, 073006 (2003), [Erratum: Phys. Rev. D73, 119901 (2006)], arXiv:hep-ph/0212229 [hep-ph] .
  * Gnendiger  _et al._ [2013] C. Gnendiger, D. Stöckinger, and H. Stöckinger-Kim, Phys. Rev. D88, 053005 (2013), arXiv:1306.5546 [hep-ph] .
  * Davier  _et al._ [2017] M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, Eur. Phys. J. C77, 827 (2017), arXiv:1706.09436 [hep-ph] .
  * Keshavarzi  _et al._ [2018] A. Keshavarzi, D. Nomura, and T. Teubner, Phys. Rev. D97, 114025 (2018), arXiv:1802.02995 [hep-ph] .
  * Colangelo  _et al._ [2019] G. Colangelo, M. Hoferichter, and P. Stoffer, JHEP 02, 006 (2019), arXiv:1810.00007 [hep-ph] .
  * Hoferichter  _et al._ [2019] M. Hoferichter, B.-L. Hoid, and B. Kubis, JHEP 08, 137 (2019), arXiv:1907.01556 [hep-ph] .
  * Davier  _et al._ [2020] M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, Eur. Phys. J. C80, 241 (2020), [Erratum: Eur. Phys. J. C80, 410 (2020)], arXiv:1908.00921 [hep-ph] .
  * Keshavarzi  _et al._ [2020] A. Keshavarzi, D. Nomura, and T. Teubner, Phys. Rev. D101, 014029 (2020), arXiv:1911.00367 [hep-ph] .
  * Kurz  _et al._ [2014] A. Kurz, T. Liu, P. Marquard, and M. Steinhauser, Phys. Lett. B734, 144 (2014), arXiv:1403.6400 [hep-ph] .
  * Melnikov and Vainshtein [2004] K. Melnikov and A. Vainshtein, Phys. Rev. D70, 113006 (2004), arXiv:hep-ph/0312226 [hep-ph] .
  * Masjuan and Sanchez-Puertas [2017] P. Masjuan and P. Sanchez-Puertas, Phys. Rev. D 95, 054026 (2017), arXiv:1701.05829 [hep-ph] .
  * Colangelo  _et al._ [2017] G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer, JHEP 04, 161 (2017), arXiv:1702.07347 [hep-ph] .
  * Hoferichter  _et al._ [2018] M. Hoferichter, B.-L. Hoid, B. Kubis, S. Leupold, and S. P. Schneider, JHEP 10, 141 (2018), arXiv:1808.04823 [hep-ph] .
  * Gérardin  _et al._ [2019] A. Gérardin, H. B. Meyer, and A. Nyffeler, Phys. Rev. D 100, 034520 (2019), arXiv:1903.09471 [hep-lat] .
  * Bijnens  _et al._ [2019] J. Bijnens, N. Hermansson-Truedsson, and A. Rodríguez-Sánchez, Phys. Lett. B798, 134994 (2019), arXiv:1908.03331 [hep-ph] .
  * Colangelo  _et al._ [2020] G. Colangelo, F. Hagelstein, M. Hoferichter, L. Laub, and P. Stoffer, JHEP 03, 101 (2020), arXiv:1910.13432 [hep-ph] .
  * Blum  _et al._ [2020] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, and C. Lehner, Phys. Rev. Lett. 124, 132002 (2020), arXiv:1911.08123 [hep-lat] .
  * Colangelo  _et al._ [2014] G. Colangelo, M. Hoferichter, A. Nyffeler, M. Passera, and P. Stoffer, Phys. Lett. B735, 90 (2014), arXiv:1403.7512 [hep-ph] .
  * Prades  _et al._ [2009] J. Prades, E. de Rafael, and A. Vainshtein, Adv. Ser. Direct. High Energy Phys. 20, 303 (2009), arXiv:0901.0306 [hep-ph] .
  * Nyffeler [2009] A. Nyffeler, Phys. Rev. D 79, 073012 (2009), arXiv:0901.1172 [hep-ph] .
  * Jegerlehner and Nyffeler [2009] F. Jegerlehner and A. Nyffeler, Phys. Rept. 477, 1 (2009), arXiv:0902.3360 [hep-ph] .
  * Jegerlehner [2018] F. Jegerlehner, EPJ Web Conf. 166, 00022 (2018), arXiv:1705.00263 [hep-ph] .
  * Green  _et al._ [2015] J. Green, O. Gryniuk, G. von Hippel, H. B. Meyer, and V. Pascalutsa, Phys. Rev. Lett. 115, 222003 (2015), arXiv:1507.01577 [hep-lat] .
  * Gérardin  _et al._ [2018] A. Gérardin, J. Green, O. Gryniuk, G. von Hippel, H. B. Meyer, V. Pascalutsa, and H. Wittig, Phys. Rev. D 98, 074501 (2018), arXiv:1712.00421 [hep-lat] .
  * Pauk and Vanderhaeghen [2014] V. Pauk and M. Vanderhaeghen, Eur. Phys. J. C 74, 3008 (2014), arXiv:1401.0832 [hep-ph] .
  * Danilkin and Vanderhaeghen [2017] I. Danilkin and M. Vanderhaeghen, Phys. Rev. D 95, 014019 (2017), arXiv:1611.04646 [hep-ph] .
  * Jegerlehner [2017] F. Jegerlehner,  _The Anomalous Magnetic Moment of the Muon_, Vol. 274 (Springer, Cham, 2017).
  * Knecht  _et al._ [2018] M. Knecht, S. Narison, A. Rabemananjara, and D. Rabetiarivony, Phys. Lett. B 787, 111 (2018), arXiv:1808.03848 [hep-ph] .
  * Eichmann  _et al._ [2020] G. Eichmann, C. S. Fischer, and R. Williams, Phys. Rev. D 101, 054015 (2020), arXiv:1910.06795 [hep-ph] .
  * Roig and Sanchez-Puertas [2020] P. Roig and P. Sanchez-Puertas, Phys. Rev. D 101, 074019 (2020), arXiv:1910.02881 [hep-ph] .
  * Gérardin  _et al._ [2016] A. Gérardin, H. B. Meyer, and A. Nyffeler, Phys. Rev. D 94, 074507 (2016), arXiv:1607.08174 [hep-lat] .
  * Alexandrou  _et al._ [2022] C. Alexandrou _et al._ , (2022), arXiv:2212.06704 [hep-lat] .
  * Burri  _et al._ [2022] S. Burri _et al._ , (2022), arXiv:2212.10300 [hep-lat] .
  * Blum  _et al._ [2014] T. Blum, S. Chowdhury, M. Hayakawa, and T. Izubuchi, (2014), arXiv:1407.2923 [hep-lat] .
  * Blum  _et al._ [2016a] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, and C. Lehner, Phys. Rev. D 93, 014503 (2016a), arXiv:1510.07100 [hep-lat] .
  * Blum  _et al._ [2017a] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, and C. Lehner, Phys. Rev. Lett. 118, 022005 (2017a), arXiv:1610.04603 [hep-lat] .
  * Hayakawa and Uno [2008] M. Hayakawa and S. Uno, Prog. Theor. Phys. 120, 413 (2008), arXiv:0804.2044 [hep-ph] .
  * Asmussen  _et al._ [2016] N. Asmussen, J. Green, H. B. Meyer, and A. Nyffeler, PoS LATTICE2016, 164 (2016), arXiv:1609.08454 [hep-lat] .
  * Blum  _et al._ [2017b] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, and C. Lehner, Phys. Rev. D 96, 034515 (2017b), arXiv:1705.01067 [hep-lat] .
  * Chao  _et al._ [2020] E.-H. Chao, A. Gérardin, J. R. Green, R. J. Hudspith, and H. B. Meyer, Eur. Phys. J. C 80, 869 (2020), arXiv:2006.16224 [hep-lat] .
  * Chao  _et al._ [2021] E.-H. Chao, R. J. Hudspith, A. Gérardin, J. R. Green, H. B. Meyer, and K. Ottnad, Eur. Phys. J. C 81, 651 (2021), arXiv:2104.02632 [hep-lat] .
  * Chao  _et al._ [2022] E.-H. Chao, R. J. Hudspith, A. Gérardin, J. R. Green, and H. B. Meyer, Eur. Phys. J. C 82, 664 (2022), arXiv:2204.08844 [hep-lat] .
  * Laporta and Remiddi [1993] S. Laporta and E. Remiddi, Phys. Lett. B 301, 440 (1993).
  * Laporta and Remiddi [1991] S. Laporta and E. Remiddi, Phys. Lett. B 265, 182 (1991).
  * Blum  _et al._ [2016b] T. Blum _et al._ (RBC, UKQCD), Phys. Rev. D 93, 074505 (2016b), arXiv:1411.7017 [hep-lat] .
  * Renfrew  _et al._ [2008] D. Renfrew, T. Blum, N. Christ, R. Mawhinney, and P. Vranas, PoS LATTICE2008, 048 (2008), arXiv:0902.2587 [hep-lat] .
  * Abramczyk  _et al._ [2017] M. Abramczyk, S. Aoki, T. Blum, T. Izubuchi, H. Ohki, and S. Syritsyn, Phys. Rev. D 96, 014501 (2017).
  * Blum  _et al._ [2018] T. Blum, P. A. Boyle, V. Gülpers, T. Izubuchi, L. Jin, C. Jung, A. Jüttner, C. Lehner, A. Portelli, and J. T. Tsang (RBC, UKQCD), Phys. Rev. Lett. 121, 022003 (2018), arXiv:1801.07224 [hep-lat] .
  * Blum  _et al._ [2023] T. Blum _et al._ , (2023), arXiv:2301.08696 [hep-lat] .
  * Clark  _et al._ [2017] M. A. Clark, C. Jung, and C. Lehner, in  _35th International Symposium on Lattice Field Theory (Lattice 2017) Granada, Spain, June 18-24, 2017_ (2017) arXiv:1710.06884 [hep-lat] .
  * Blum  _et al._ [2013] T. Blum, T. Izubuchi, and E. Shintani, Phys. Rev. D 88, 094503 (2013), arXiv:1208.4349 [hep-lat] .
  * Detmold  _et al._ [2021] W. Detmold, D. J. Murphy, A. V. Pochinsky, M. J. Savage, P. E. Shanahan, and M. L. Wagman, Phys. Rev. D 104, 034502 (2021), arXiv:1908.07050 [hep-lat] .
  * Li  _et al._ [2021] Y. Li, S.-C. Xia, X. Feng, L.-C. Jin, and C. Liu, Phys. Rev. D 103, 014514 (2021), arXiv:2009.01029 [hep-lat] .
  * Bijnens and Relefors [2016] J. Bijnens and J. Relefors, JHEP 09, 113 (2016), arXiv:1608.01454 [hep-ph] .
  * Jin  _et al._ [2016] L. Jin, T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, C. Jung, and C. Lehner, PoS LATTICE2016, 181 (2016), arXiv:1611.08685 [hep-lat] .
  * Note [1] One can use different setup for the connected and disconnected diagrams, in which case this relation will not hold exactly. For example, Ref. [48] uses different setups.
  * Zhao [2022] Y. Zhao,  _Lattice Calculation of the $\pi^{0}\rightarrow e^{+}e^{-}$ and the $K_{L}\rightarrow\gamma\gamma$ Decays_, Ph.D. thesis, Columbia U. (2022).
  * Asmussen  _et al._ [2022] N. Asmussen, E.-H. Chao, A. Gérardin, J. R. Green, R. J. Hudspith, H. B. Meyer, and A. Nyffeler, (2022), arXiv:2210.12263 [hep-lat] .
  * Feng  _et al._ [2012] X. Feng, S. Aoki, H. Fukaya, S. Hashimoto, T. Kaneko, J.-i. Noaki, and E. Shintani, Phys. Rev. Lett. 109, 182001 (2012), arXiv:1206.1375 [hep-lat] .
  * Moussallam [1995] B. Moussallam, Phys. Rev. D 51, 4939 (1995), arXiv:hep-ph/9407402 .
  * Knecht  _et al._ [1999] M. Knecht, S. Peris, M. Perrottet, and E. de Rafael, Phys. Rev. Lett. 83, 5230 (1999), arXiv:hep-ph/9908283 .
  * Ramsey-Musolf and Wise [2002] M. J. Ramsey-Musolf and M. B. Wise, Phys. Rev. Lett. 89, 041601 (2002), arXiv:hep-ph/0201297 .
  * Kinoshita  _et al._ [1985] T. Kinoshita, B. Nizic, and Y. Okamoto, Phys. Rev. D 31, 2108 (1985).
  * Hayakawa  _et al._ [1995] M. Hayakawa, T. Kinoshita, and A. I. Sanda, Phys. Rev. Lett. 75, 790 (1995), arXiv:hep-ph/9503463 .
  * Bijnens [2016] J. Bijnens, EPJ Web Conf. 118, 01002 (2016), arXiv:1510.05796 [hep-ph] .
  * Christ  _et al._ [2022] N. Christ, X. Feng, L. Jin, C. Tu, and Y. Zhao, (2022), arXiv:2208.03834 [hep-lat] .