# Lattice calculation of the $\pi^{0}$, $\eta$ and $\eta^{\prime}$ transition
form factors and the hadronic light-by-light contribution to the muon $g-2$

Antoine Gérardin<EMAIL_ADDRESS>Aix-Marseille Université,
Université de Toulon, CNRS, CPT, Marseille, France  Willem E. A. Verplanke
<EMAIL_ADDRESS>Aix-Marseille Université, Université de
Toulon, CNRS, CPT, Marseille, France  Gen Wang  Aix-Marseille Université,
Université de Toulon, CNRS, CPT, Marseille, France  Zoltan Fodor  Department
of Physics, University of Wuppertal, D-42119 Wuppertal, Germany  Jülich
Supercomputing Centre, Forschungszentrum Jülich, D-52428 Jülich, Germany
Institute for Theoretical Physics, Eötvös University, H-1117 Budapest, Hungary
Department of Physics, Pennsylvania State University, University Park, PA
16802, USA  Jana N. Guenther  Department of Physics, University of Wuppertal,
D-42119 Wuppertal, Germany  Laurent Lellouch  Aix-Marseille Université,
Université de Toulon, CNRS, CPT, Marseille, France  Kalman K. Szabo
Department of Physics, University of Wuppertal, D-42119 Wuppertal, Germany
Jülich Supercomputing Centre, Forschungszentrum Jülich, D-52428 Jülich,
Germany  Lukas Varnhorst  Department of Physics, University of Wuppertal,
D-42119 Wuppertal, Germany

###### Abstract

In this paper we present a first ab-initio calculation of the $\pi^{0}$,
$\eta$ and $\eta^{\prime}$ transition form factors performed with physical
light-quark masses. We provide a complete parametrization of the form factors
that includes both single and double-virtual kinematics. Our results are
compared with experimental measurements of the form factors in the space-like
region and with the measured two-photon decay widths. In a second step, our
parametrizations of the transition form factors are used to compute the
dominant pseudoscalar-pole contributions to the hadronic light-by-light
scattering in the muon $g-2$. Our final result reads $a_{\mu}^{\rm hlbl,ps-
pole}=(85.1\pm 5.2)\times 10^{-11}$. Although the pion-pole is dominant, we
confirm that, together, the $\eta$ and $\eta^{\prime}$ provide roughly half of
its contribution.

##  I Introduction

Figure 1: Pseudoscalar-pole contribution to the hadronic light-by-light
scattering diagram. The blob on the left of the equality represents the full
HLbL four-point function. The blobs on the right-hand side represent the
transition form factors of the light pseudoscalar mesons. The solid, dashed
and wavy lines represent the muon, the intermediate pseudoscalar mesons, and
the photons, respectively.

The combined experimental result for the anomalous magnetic moment of the
muon, from the Fermilab collaboration [1] and the E821 experiment at
Brookhaven [2], has reached a precision of 0.35 ppm. It translates to a
4.2$\sigma$ tension with the theory average, as published by the $g-2$ theory
initiative [3]. The theory error is completely dominated by effects of the
strong interaction which can be separated into two distinct contributions: the
hadronic vacuum polarization (HVP) which enters at order $\alpha^{2}$ and the
hadronic light-by-light scattering (HLbL) at order $\alpha^{3}$. The
experimental uncertainties [4, 5] are expected to be reduced by a further
factor of 3 by the Fermilab experiment [4] and a new experiment at J-PARC [5]
is under construction. To fully benefit from these two new experiments, a
precision of 0.2% for the HVP and below 10% for the HLbL is required.

Over the past few years, lattice QCD has made significant progress toward
reducing the uncertainties on the HVP contribution [6, 7, 8, 9, 10, 11, 12,
13, 14, 15, 16, 17, 18, 19, 20] and the Budapest-Marseille-Wuppertal
collaboration published the first sub-percent calculation [16]. Interestingly,
this result is in slight tension, at the level of 2.1$\sigma$, with the data-
driven estimate [21, 22, 23, 24, 25, 26] and tends to reduce the discrepancy
between the theory estimate and the world-average measurement to the level of
1.5$\sigma$. Significant attention has been devoted to clarify this
discrepancy [27, 28, 29, 30]. The situation for the HLbL contribution is
different. Two groups have presented direct lattice calculations where the
hadronic correlation function of four electromagnetic currents is computed on
the lattice [31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42]. This correlator,
computed in position-space, is weighted by a kernel function that represents
the QED part of the diagram (see Fig. 1). Beside the use of different fermion
formulations, the main difference between the two groups lies in the treatment
of the weight function that differ between the Mainz group [41] and the
RBC/UKQCD collaboration [35]. Both lattice estimates are compatible with each
other and with the data-driven result based on the dispersive framework
developed in [43, 44, 45, 46, 47]. However, the current precision of this
comparison is at the level of 15% and more stringent tests require improvement
in both frameworks.

On the lattice, the main challenge consists in reducing the statistical error
that increases rapidly at long distances and to control possibly large finite-
volume corrections. Fortunately, both sources of uncertainty can be strongly
suppressed if one knows the pion transition form factor (TFF) on the same set
of gauge ensembles [48]. The knowledge of the $\eta$ and $\eta^{\prime}$ TFFs
would allow to test the saturation of the signal at even shorter distances.

An alternative to the direct lattice calculation, is the data driven approach
[43, 44, 45, 46, 47]. It is based on a dispersive framework, similar to the
one used for the HVP, but is much more complicated due to the complex
analytical structure of the four-point function involved in the HLbL diagram.
It is expected that, among all possible intermediate states, the largest
contributions are given by a handful of states: the light pseudoscalar mesons
($\pi^{0}$, $\eta$ and $\eta^{\prime}$). In this dispersive framework, the
contribution from the pseudoscalar states can be unambiguously computed once
the TFFs, describing the interaction of the mesons with two virtual photons,
are known (see Fig. 1).

Although the form factors are, in principle, accessible through experiments,
their measurement in the whole kinematic range relevant for the $g-2$ is
challenging. Only the single-virtual form factor has been measured by the
BaBar [49], CELLO [50], CLEO [51] and Belle [52] collaborations for space-like
virtualities above 0.6 $\mathrm{GeV}^{2}$. Recently, the BES III collaboration
has presented preliminary results down to $0.3\leavevmode\nobreak\
\mathrm{GeV}^{2}$. The single-virtual $\eta$ and $\eta^{\prime}$ form factors
have been measured in the space-like region by CELLO [50], CLEO [51], BaBar
[53] and the L3 experiment at LEP [54], covering the region from
$0.3\leavevmode\nobreak\ \mathrm{GeV}^{2}$ up to $40\leavevmode\nobreak\
\mathrm{GeV}^{2}$. The BaBar collaboration has published first results for the
double-virtual TFF associated with the $\eta^{\prime}$ meson, but only for
virtualities larger than $Q^{2}\approx 6.5\leavevmode\nobreak\
\mathrm{GeV}^{2}$ [55]. When both photons are real, the form factors are
related to the two-photon decay width. The latter has been measured for all
three mesons [56, 57] and is an important constraint to the low-energy
behavior of the TFFs.

In [58, 59] it has been shown that a precise lattice calculation of the pion
TFF in the kinematical range relevant for the $g-2$ is possible, leading to
$a_{\mu}^{\mathrm{HLbL};\pi^{0}}=(59.7\pm 3.6)\times 10^{-11}$. This result is
in good agreement with the data-driven estimate [60, 61]. Using similar
techniques, the $\eta$ pole contribution has been computed by the ETM
collaboration, but at a single lattice spacing [62]. In this work, we provide
a new estimate of the pion-pole contribution, using a different fermionic
action as compared to [59], and we extend the calculation to the $\eta$ and
$\eta^{\prime}$ mesons. We provide parametrizations of the form factors at the
physical point that are then used to estimate the light pseudoscalar-pole
contributions to the HLbL diagram. To this end, we are preparing a dedicated
spectroscopy analysis of the $\eta$ and $\eta^{\prime}$ mesons using staggered
quarks on the same set of ensembles [63].

This paper is organized as follows. In Section II we describe our strategy to
extract the TFFs in Euclidean spacetime using lattice QCD simulations. In
Section III we present our lattice setup and explain how the relevant
correlation functions are computed with staggered quarks. Then, in Section IV,
we present our analysis and we discuss the main sources of systematic
uncertainties. Finally, in Section V, we present our results for all three
TFFs and we discuss phenomenological applications. Tables with the
coefficients of the TFF parametrizations in the continuum limit, and the
associated covariance matrices, are provided in Appendix A. In Appendix B, we
compare different estimators to compute single vector current insertions in a
quark loop.

##  II Methodology

###  II.1 Form factors

In this section, we follow the notations introduced in [58, 59]. The TFF
describing the interaction between a pseudoscalar meson
$P=\pi^{0},\eta,\eta^{\prime}$, with momentum $p$, and two off-shell photons,
with momenta $q_{1}$ and $q_{2}$, such that $p=q_{1}+q_{2}$, is defined via
the following matrix element

$M_{\mu\nu}(p,q_{1})=i\int\mathrm{d}^{4}x\,e^{iq_{1}\cdot
x}\,\langle\Omega|T\\{J_{\mu}(x)J_{\nu}(0)\\}|P(p)\rangle=\epsilon_{\mu\nu\alpha\beta}\,q_{1}^{\alpha}\,q_{2}^{\beta}\,{\cal
F}_{P\gamma^{*}\gamma^{*}}(q_{1}^{2},q_{2}^{2})\,.$ (1)

In this equation, $J_{\mu}$ is the hadronic component of the electromagnetic
current and $\epsilon_{\mu\nu\alpha\beta}$ is the fully antisymmetric tensor
with $\epsilon^{0123}=1$. We have implicitly treated the meson as an
asymptotic state. This is true for the pion but also for the $\eta$ meson in
the isospin limit of QCD (the decays $\eta\to\pi^{+}\pi^{-}\pi^{0}$ and
$\eta\to 3\pi^{0}$ break isospin symmetry). On the other hand, the
$\eta^{\prime}$ meson can decay via the strong interaction, even in the
isospin limit, and its main decay modes are
$\eta^{\prime}\to\eta\pi^{+}\pi^{-}$ and $\eta^{\prime}\to\eta\pi^{0}\pi^{0}$.
However, the smallness of the decay width justifies the narrow-width
approximation in which the meson is treated as stable particle. Below the
threshold for hadron production in the vector channel111In the isovector case,
the treshold $s_{0}$ is given by $4m_{\pi}^{2}$., the TFF is obtained in
Euclidean spacetime after analytical continuation [64, 65]. Using the
superscript E, the Euclidean matrix element reads

$M_{\mu\nu}=(i^{n_{0}})M_{\mu\nu}^{\rm E},\quad M_{\mu\nu}^{\rm
E}\equiv-\int_{-\infty}^{+\infty}\mathrm{d}\tau\,e^{\omega_{1}\tau}\int\mathrm{d}^{3}z\,e^{-i\vec{q}_{1}\cdot\vec{z}}\,\langle
0|T\left\\{J_{\mu}(\vec{z},\tau)J_{\nu}(\vec{0},0)\right\\}|P(p)\rangle\,,$
(2)

where $n_{0}$ denotes the number of temporal indices carried by the two vector
currents. Different virtualities can be reached by tuning the real, free
parameter $\omega_{1}$ such that $q_{1}=(\omega_{1},\vec{q}_{1})$ and
$q_{2}=p-q_{1}$. It is convenient to write Eq. (2) as an integral over $\tau$,
the time separation between the two vector currents222Note that our
conventions slightly differ from [59].

$M_{\mu\nu}^{\rm
E}=\int_{-\infty}^{\infty}\,\mathrm{d}\tau\,\widetilde{A}^{(P)}_{\mu\nu}(\tau)\,e^{\omega_{1}\tau}\,.$
(3)

Since we are working in the isospin limit with equal-mass $u$ and $d$ quarks,
the pion does not mix with the $\eta$ and $\eta^{\prime}$ mesons and an
interpolating operator for this meson is given by the pseudoscalar density
$P_{3}(x)$ below. If we further assume SU(3) flavor symmetry then the $\eta$
and $\eta^{\prime}$ mesons can be identified with the octet $\eta_{8}$ and
singlet $\eta_{0}$ pseudoscalar mesons. Those states can be studied using the
densities $P_{8}(x)$ and $P_{0}(x)$ respectively. We have:

$\displaystyle P_{3}(x)$
$\displaystyle=\frac{1}{\sqrt{2}}\left(\overline{u}\gamma_{5}u(x)-\overline{d}\gamma_{5}d(x)\right),$
(4a) $\displaystyle P_{8}(x)$
$\displaystyle=\frac{1}{\sqrt{6}}\left(\overline{u}\gamma_{5}u(x)+\overline{d}\gamma_{5}d(x)-2\overline{s}\gamma_{5}s(x)\right),$
(4b) $\displaystyle P_{0}(x)$
$\displaystyle=\frac{1}{\sqrt{3}}\left(\overline{u}\gamma_{5}u(x)+\overline{d}\gamma_{5}d(x)+\overline{s}\gamma_{5}s(x)\right).$
(4c)

Away from the SU(3) flavor limit, the $\eta_{8}$ and $\eta_{0}$ states mix to
give the physical $\eta$ and $\eta^{\prime}$ mesons that have non-vanishing
overlaps with both $P_{8}$ and $P_{0}$ . We now explain how this effect is
taken into account in this work.

To extract the matrix element in Eq. (2), we consider the following three-
point correlation functions

$C^{(i)}_{\mu\nu}(\tau,t_{P})=\sum_{\vec{x},\vec{z}}\langle
J_{\mu}(\vec{z},t_{i})J_{\nu}(\vec{0},t_{f})\mathcal{O}_{i}(\vec{x},t_{0})\rangle\,e^{i\vec{p}\cdot\vec{x}}e^{-i\vec{q}_{1}\cdot\vec{z}}\,,$
(5)

where $t_{P}={\rm min}(t_{f}-t_{0},t_{i}-t_{0})$ is the minimal time
separation between the pseudoscalar interpolating operator and either of the
two vector currents and $\tau=t_{i}-t_{f}$. The lattice operators
$\mathcal{O}_{i}$ have the same quantum numbers as the pseudoscalar densities
$P_{i}$. The precise definition of those operators will be given in Section
III.2. For the pion, neglecting excited states and wrap around contributions,
we obtain the asymptotic behavior

$\displaystyle C^{(3)}_{\mu\nu}(\tau,t_{P})=$
$\displaystyle\sum_{\vec{z}}\langle
0|J_{\mu}(\vec{z},\tau)J_{\nu}(\vec{0},0)|\pi(\vec{p})\rangle
e^{-i\vec{q}_{1}\cdot\vec{z}}\times\frac{1}{2E_{\pi}}\langle\pi(\vec{p})|\mathcal{O}_{3}|0\rangle
e^{E_{\pi}(t_{0}-t_{f})}$ (6)

where $E_{\pi}$ is the energy of the pion and
$Z_{\pi}=\langle\pi|\mathcal{O}_{3}|0\rangle$ is the overlap factor of the
pseudoscalar operator with the pion state. It leads to

$\widetilde{A}^{(\pi)}_{\mu\nu}(\tau)=\lim_{t_{P}\rightarrow+\infty}\widetilde{A}_{\mu\nu}^{(\pi);{\rm
eff}}(\tau,t_{P})\,,\quad\widetilde{A}_{\mu\nu}^{(\pi);{\rm
eff}}(\tau,t_{P})=\frac{2E_{\pi}}{Z_{\pi}}e^{E_{\pi}(t_{f}-t_{0})}C^{(3)}_{\mu\nu}(\tau,t_{P})\,.$
(7)

A similar expression could be used for the $\eta$ meson, looking at the
correlator $C^{(8)}_{\mu\nu}(\tau,t_{P})$. Indeed, the correlation function
$C^{(8)}_{\mu\nu}(\tau,t_{P})$ is expected to be dominated by the $\eta$
meson, with little contribution from the heavier $\eta^{\prime}$ that can be
seen as an excited state. This is the procedure followed in [62]. However,
this method has some disadvantages. First, the long time asymptotic is reached
at large values of $t_{P}$ where the signal-to-noise ratio deteriorates
rapidly. Second, this method is not applicable to the $\eta^{\prime}$ state. A
comparison of two approaches for the $\eta$ meson is presented in Section
IV.4.

Instead, the $\eta$ and $\eta^{\prime}$ TFFs are extracted from a $2\times 2$
matrix of three-point functions. Neglecting higher excited states, the
spectral decomposition for $i=0,8$ now reads

$\displaystyle C^{(i)}_{\mu\nu}(\tau,t_{P})=$
$\displaystyle\sum_{\vec{z}}\langle
0|J_{\mu}(\vec{z},\tau)J_{\nu}(\vec{0},0)|\eta(\vec{p})\rangle
e^{-i\vec{q}_{1}\cdot\vec{z}}\times\frac{1}{2E_{\eta}}\langle\eta(\vec{p})|\mathcal{O}_{i}|0\rangle
e^{E_{\eta}(t_{0}-t_{f})}$ $\displaystyle+$
$\displaystyle\sum_{\vec{z}}\langle
0|J_{\mu}(\vec{z},\tau)J_{\nu}(\vec{0},0)|\eta^{\prime}(\vec{p})\rangle
e^{-i\vec{q}_{1}\cdot\vec{z}}\times\frac{1}{2E_{\eta^{\prime}}}\langle\eta^{\prime}(\vec{p})|\mathcal{O}_{i}|0\rangle
e^{E_{\eta}^{\prime}(t_{0}-t_{f})}\,.$ (8)

Or, in a more compact matrix notation

$\begin{pmatrix}C_{\mu\nu}^{(8)}\\\\[4.2679pt]
C_{\mu\nu}^{(0)}\end{pmatrix}=\begin{pmatrix}T^{(8)}_{\eta}&T^{(8)}_{\eta^{\prime}}\\\
T^{(0)}_{\eta}&T^{(0)}_{\eta^{\prime}}\end{pmatrix}\begin{pmatrix}\widetilde{A}^{(\eta)}_{\mu\nu}\\\\[3.69885pt]
\widetilde{A}^{(\eta^{\prime})}_{\mu\nu}\end{pmatrix}\,,$ (9)

with

$T_{n}^{(i)}=\frac{Z_{n}^{(i)}}{2E_{n}}e^{-E_{n}(t_{f}-t_{0})}\,,\quad\widetilde{A}^{(n)}_{\mu\nu}(\tau)=\sum_{\vec{z}}\langle
0|J_{\mu}(\vec{z},\tau)J_{\nu}(\vec{0},0)|n(\vec{p})\rangle\,e^{-i\vec{q}_{1}\cdot\vec{z}}\,.$
(10)

The four overlap factors are defined by $Z_{n}^{(i)}=\langle
n|\mathcal{O}_{i}|0\rangle$. Inverting Eq. (9), one obtains the $\eta$ and
$\eta^{\prime}$ TFFs in terms of the correlators computed on the lattice

$\displaystyle\widetilde{A}^{(\eta)}_{\mu\nu}$
$\displaystyle=\cos^{2}\phi_{I}\
\frac{C_{\mu\nu}^{(8)}}{T_{\eta}^{(8)}}+\sin^{2}\phi_{I}\
\frac{C_{\mu\nu}^{(0)}}{T_{\eta}^{(0)}}$ (11a)
$\displaystyle\widetilde{A}^{(\eta^{\prime})}_{\mu\nu}$
$\displaystyle=\sin^{2}\phi_{I}\
\frac{C_{\mu\nu}^{(8)}}{T_{\eta^{\prime}}^{(8)}}+\cos^{2}\phi_{I}\
\frac{C_{\mu\nu}^{(0)}}{T_{\eta^{\prime}}^{(0)}}$ (11b)

where the mixing angle in the isospin basis, $\phi_{I}$, is given by
$\tan^{2}\phi_{I}=-(Z^{(8)}_{\eta^{\prime}}Z^{(0)}_{\eta})/(Z^{(8)}_{\eta}Z^{(0)}_{\eta^{\prime}})$.

In practice, the masses and overlap factors $Z_{n}^{(i)}$ can be extracted by
solving the Generalized Eigenvalue Problem (GEVP) using a matrix of
pseudoscalar two-point correlation functions

$C^{\rm 2pt}(t)v_{n}(t,t_{0})=\lambda_{n}(t,t_{0})C^{\rm
2pt}(t_{0})v_{n}(t,t_{0})\,,$ (12)

where $C^{\rm 2pt}_{ij}(t)=\langle\mathcal{O}_{i}(t)\mathcal{O}_{j}(0)\rangle$
with $i,j\in 0,8$ and with $t_{0}$ a free parameter. The effective masses can
be extracted from the logarithmic derivative of the eigenvalues $\lambda_{n}$
[66]

$am_{n}^{\rm eff}(t)=\log\frac{\lambda_{n}(t,t_{0})}{\lambda_{n}(t+a,t_{0})}$
(13)

that can be eventually fitted to a constant at large times. The overlap
factors can be obtained from the eigenvectors $v_{n}$ through

$Z_{n}^{{(i)}\,{\rm eff}}(t)=\sqrt{2E_{n}}\ \sum_{j}C^{\rm
2pt}_{ij}(t)\,v_{nj}(t,t_{0})\times\left(\frac{\lambda_{n}(t,t_{0})}{\lambda_{n}(t+a,t_{0})}\right)^{t/a-t_{0}/(2a)}$
(14)

assuming the normalization

$v_{m}(t,t_{0})^{\rm T}\,C^{\rm 2pt}(t_{0})\,v_{n}(t,t_{0})=\delta_{n,m}\,.$
(15)

However, we prefer to follow the strategy presented in [63] where a fit of the
pseudoscalar two-point correlation matrix is performed. This method leads to
numerically more stable results.

###  II.2 Pseudoscalar-pole contribution to HLbL scattering in $(g-2)_{\mu}$

Table 1: Parameters of the simulations: the bare coupling $\beta=6/g_{0}^{2}$, the lattice dimensions $L^{3}\times T$, the lattice spacing $a$ and the spatial extent $L$ in physical units, the bare light and strange quark masses, the number of gauge configurations. The last three columns indicate which TFFs have been computed. $\quad\beta\quad$ | $L^{3}\times T$ | $a\leavevmode\nobreak\ [\mathrm{fm}]$ |  $L\leavevmode\nobreak\ [\mathrm{fm}]$ | $am_{l}$ | $am_{s}$ |  $\\#$ cnfgs | $\pi^{0}$ | $\eta$ | $\eta^{\prime}$
---|---|---|---|---|---|---|---|---|---
$3.7000$ | $24^{3}\times 48$ | 0.1315 | 3.2 | 0.00205349 | 0.0572911 | 700 | ✓ | ✓ | ✓
| $32^{3}\times 64$ |  | 4.2 | 0.00205349 | 0.0572911 | 900 | ✓ | ✓ | ✓
| $48^{3}\times 64$ |  | 6.3 | 0.00205349 | 0.0572911 | 900 | ✓ | ✓ | ✓
$3.7500$ | $56^{3}\times 96$ | 0.1191 | 6.7 | 0.00184096 | 0.0495930 | 500 | ✓ | ✓ | $-$
| $56^{3}\times 96$ |  | 6.7 | 0.00176877 | 0.0516173 | 500 | ✓ | ✓ | $-$
| $56^{3}\times 96$ |  | 6.7 | 0.00184096 | 0.0516173 | 500 | ✓ | ✓ | $-$
$3.7753$ | $28^{3}\times 56$ | 0.1116 | 3.1 | 0.00171008 | 0.0476146 | 850 | ✓ | ✓ | ✓
| $56^{3}\times 84$ |  | 6.4 | 0.00171008 | 0.0476146 | 500 | ✓ | ✓ | $-$
| $56^{3}\times 84$ |  | 6.4 | 0.00171008 | 0.0485669 | 500 | ✓ | ✓ | $-$
| $56^{3}\times 84$ |  | 6.4 | 0.00174428 | 0.0461862 | 500 | ✓ | ✓ | ✓
$3.8400$ | $32^{3}\times 64$ | 0.0952 | 3.0 | 0.00151556 | 0.0431935 | 1100 | ✓ | ✓ | ✓
| $32^{3}\times 64$ |  | 3.0 | 0.00143 | 0.0431935 | 1050 | ✓ | ✓ | ✓
| $32^{3}\times 64$ |  | 3.0 | 0.001455 | 0.04075 | 1100 | ✓ | ✓ | ✓
| $32^{3}\times 64$ |  | 3.0 | 0.001455 | 0.03913 | 1100 | ✓ | ✓ | ✓
| $64^{3}\times 96$ |  | 6.1 | 0.00151556 | 0.0431935 | 500 | ✓ | $-$ | $-$
| $64^{3}\times 96$ |  | 6.1 | 0.001455 | 0.04075 | 1100 | ✓ | $-$ | $-$
$3.9200$ | $40^{3}\times 80$ | 0.0787 | 3.1 | 0.001207 | 0.032 | 550 | ✓ | ✓ | ✓
| $40^{3}\times 80$ |  | 3.1 | 0.0012 | 0.0332856 | 350 | ✓ | ✓ | ✓
| $80^{3}\times 128$ |  | 6.3 | 0.001172 | 0.03244 | 500 | ✓ | $-$ | $-$
| $80^{3}\times 128$ |  | 6.3 | 0.0012 | 0.0332856 | 500 | ✓ | $-$ | $-$
$4.0126$ | $48^{3}\times 96$ | 0.0640 | 3.1 | 0.00095897 | 0.0264999 | 850 | ✓ | ✓ | ✓
| $48^{3}\times 96$ |  | 3.1 | 0.001002 | 0.027318 | 450 | ✓ | ✓ | ✓
| $96^{3}\times 144$ |  | 6.1 | 0.000977 | 0.0264999 | 500 | ✓ | $-$ | $-$
| $96^{3}\times 144$ |  | 6.1 | 0.001002 | 0.027318 | 450 | ✓ | $-$ | $-$

The TFFs are the key ingredients to compute the pseudoscalar-pole contribution
to the HLbL diagram in the dispersive framework introduced in [43, 44, 45, 46,
47]. The master equation, first derived in Ref. [67], holds exactly for the
pion and the $\eta$ in the isospin limit and for the $\eta^{\prime}$ in the
narrow-width approximation. Starting from the two-loop integral (see Fig. 1),
all angular integrals, except one, are performed using the Gegenbauer
polynomial technique leading to

$a_{\mu}^{\mathrm{HLbL};P}=\left(\frac{\alpha}{\pi}\right)^{3}\int_{0}^{\infty}\\!\\!\\!dQ_{1}\\!\\!\int_{0}^{\infty}\\!\\!\\!dQ_{2}\\!\\!\int_{-1}^{1}\\!\\!d\beta\,\Big{(}w_{1}(Q_{1},Q_{2},\beta)\,{\cal
F}_{P\gamma^{*}\gamma^{*}}(-Q_{1}^{2},-(Q_{1}+Q_{2})^{2})\times\\\ {\cal
F}_{P\gamma^{*}\gamma^{*}}(-Q_{2}^{2},0)+w_{2}(Q_{1},Q_{2},\beta)\,{\cal
F}_{P\gamma^{*}\gamma^{*}}(-Q_{1}^{2},-Q_{2}^{2})\,{\cal
F}_{P\gamma^{*}\gamma^{*}}(-(Q_{1}+Q_{2})^{2},0)\Big{)}\,,$ (16)

where $w_{1}$ and $w_{2}$ are two model-independent weight functions, $Q_{1}$
and $Q_{2}$ are the norms of the two Euclidean four-momentum vectors and
$\beta=\cos\theta$ is related to the angle between them, $Q_{1}\cdot
Q_{2}=Q_{1}Q_{2}\cos\theta$. The weight functions and the shape of the
integrand will be further discussed in Section V.4, along with the
presentation of the lattice data.

##  III Simulation details

###  III.1 Lattice ensembles

This work is based on a subset of gauge ensembles generated by the Budapest-
Marseille-Wuppertal collaboration [16]. They have been generated using
$N_{f}=2+1+1$ dynamical staggered fermions with four steps of stout smearing.
The bare quark masses have been tuned such that the Goldstone mesons are at
nearly physical pion and kaon mass. The lattice spacing is set using the
$\Omega$ baryon mass and we exploit six values of the lattice spacing in the
range $[0.0640-0.1315]$ fm to extrapolate our result to the continuum limit.
We also consider $L=3$, $4$ and $6$ fm boxes for finite-size effect studies.
Simulations are performed in the isospin limit where $m_{u}=m_{d}\equiv
m_{\ell}$. More details about those ensembles can be found in [16] and the
main properties relevant to this work are summarized in Table 1.

###  III.2 Staggered correlation functions

The extraction of the TFFs requires the calculation of the three-point
functions in Eq. (5) and the corresponding pseudoscalar two-point functions to
extract the pseudoscalar masses and overlap factors that appear in Eq. (11).
We now describe the construction of those correlators.

In the staggered formalism, the operators are chosen such that the overall
correlation function is a taste singlet. Since we are interested in the
extraction of the $\eta^{(\prime)}$ TFFs, both the two- and three-point
correlation functions involve Wick contractions with quark-disconnected
contributions for the pseudoscalar density. Such contributions require taste-
singlet pseudoscalar operators. Using the notations of [68], where operators
are given in the staggered spin-taste basis as $\Gamma_{S}\otimes\Gamma_{T}$,
the operator would correspond to a $\gamma_{5}\otimes 1$ operator. In [63] we
consider two operators [69]: a time-local 3-link operator and a time non-local
4-link operator where the operator has support on two timeslices. In this
work, we use the 4-link operator that has the advantage to highly suppress the
parity partner state such that no oscillation in time is observed in the
analysis of the two-point correlation functions. It makes the subsequent
analysis of the two and three point functions much simpler.

For the vector current, we also use a taste singlet operator. More precisely,
we use the conserved one-link operators

$J_{\mu}(x)=-\frac{1}{2}\eta_{\mu}(x)\left[\overline{\chi}(x+a\hat{\mu})U_{\mu}^{{\dagger}}(x)\chi(x)+\overline{\chi}(x)U_{\mu}(x)\chi(x+a\hat{\mu})\right]\,,$
(17)

where $\chi$ is a single component staggered fermion field and $\eta_{\mu}$ is
a phase factor (we use the convention of [68], also used in [63]). The gauge
links $U_{\mu}(x)$ ensure that correlation functions are gauge invariant. In
Eq. (17), the current is defined for a single flavor and the electromagnetic
charge factors are not included. Our choice is motivated by the presence of
disconnected contributions that involve a single vector loop and also require
a taste-singlet operator. Moreover, this current does not require
multiplicative renormalization. We note that this operator was already used in
[16] to compute the LO-HVP contribution.

Among the 16 tastes, the taste-singlet pion is the heaviest. Taste splitting
is expected to decrease as O($a^{2}$) as we approach the continuum limit [70,
71]. Between our coarsest and finest lattice spacings the taste-singlet pion
mass decreases from about 400 MeV to 200 MeV while the goldstone pion
$\gamma_{5}\otimes\gamma_{5}$ is kept at the physical pion mass: a strong
dependence on the lattice spacing is possible.

###  III.3 Computational strategy

Figure 2: The four different Wick contraction topologies that are present in
the three-point correlation function. Some topologies contains several
diagrams that are not shown for brevity. The blue and green blobs represent a
pseudoscalar or a vector current insertion respectively. Only the first two
contractions contribute to the pion TFF.

The three-point correlation function in Eq. (5) receives potential
contributions from four Wick contractions, shown in Fig. 2. In this section,
the pseudoscalar and vector currents are defined for a single flavor and the
correct charges factors associated with the operators in Eq. (4) need to be
included when building the $\pi^{0}$, $\eta$ and $\eta^{\prime}$ correlators.

The connected contribution to the three-point function is computed using
similar techniques to the one used in [58, 59]. First, a point source is
created at a random position on the timslice $t_{f}$. For each value of the
meson momentum $\vec{p}$, the pseudoscalar operator is applied to the solution
vector, restricted to the timeslice $t_{0}$, and then used as a source for a
sequential inversion. The correlation function is finally obtained by
contracting the original and the sequential solution vectors with appropriate
phase factors and shifts on timeslice $t_{i}$. This method allows us to access
many photon virtualities $\vec{q}_{1}$ and $\vec{q}_{2}$ at no additional
cost. Only additional pion momenta, or additional values of $t_{P}$, defined
after Eq. (5), require new evaluations of the quark propagator.

The main computational challenge is to obtain a good signal for the quark-
disconnected contributions. First, we need to compute the pseudoscalar and
vector loop functions

$\displaystyle L_{P}^{(f)}(x_{0},\vec{p})$
$\displaystyle=-\left(\frac{a}{L}\right)^{3}\sum_{\vec{x}}\mathrm{Tr}\left[\mathcal{O}_{P}\,S_{f}(x,x)\,\right]\,e^{i\vec{p}\cdot\vec{x}}$
(18a) $\displaystyle L_{V;\mu}^{(f)}(x_{0};\vec{q}_{1})$
$\displaystyle=-\frac{1}{2}\left(\frac{a}{L}\right)^{3}\sum_{\vec{x}}\,\eta_{\mu}(x)\,\mathrm{Tr}\left[S_{f}(x,x+a\hat{\mu})U_{\mu}^{{\dagger}}(x)+S_{f}(x+a\hat{\mu},x)U_{\mu}(x)\right]\,e^{i\vec{q}_{1}\cdot\vec{x}}$
(18b)

with $S_{f}$ the quark propagator of flavor $f=(l,s)$, $\mathcal{O}_{P}$ the
pseudoscalar taste-singlet 4-link operator presented in [63] and where the
trace runs over color indices. Those functions have been stored on disk for
all time slices and all required values of the momenta. The pseudoscalar loops
have been computed using a combination of low-mode averaging (LMA) [72, 73]
and all-mode-averaging (AMA) [74, 75, 76] combined with the Venkataraman-
Kilcup variance reduction (VKVR) trick [77] to reach the intrinsic gauge
noise. More details are provided in Ref. [63]. For the vector loops, the VKVR
trick is not applicable. Instead, since we are only interested in the light
minus strange combination $L_{V}^{(l-s)}$ (this is the only combination that
appears in the isospin limit of QCD), we adopt the split-even stochastic
estimator introduced in [78] combined with LMA and AMA to reach the gauge
noise (see Appendix B for more details).

Since we want to access many photon virtualities, the vector-vector two-point
correlation functions with flavor $f$, and denoted
$D^{(f)}_{VV;\mu\nu}(t_{f},t_{i};\vec{p},\vec{q}_{1})$ (third diagram in Fig.
2), has been computed using point sources for all source ($t_{i}$) and sink
($t_{f}$) time slices. It is a noisy contribution and we found that it was
optimal to use AMA with 64 sloppy inversions per exact solve for the light
quark and 16 for the strange quark. Sloppy inversions are defined by a fixed
number of iterations, chosen such that the bias correction is small as
compared to the statistical error. In practice, we used 1, 2 or 3 exact solves
per time slice depending on the ensemble. The corresponding three-point
correlators, with all flavor combinations, are finally obtained through

$C^{(3)}_{\mathrm{P}_{f_{2}}-\mathrm{VV}_{f_{1}};\mu\nu}(\tau,t_{P};\vec{p},\vec{q}_{1})=\frac{a}{T}\sum_{t_{i}}\
\langle
D^{(f_{1})}_{VV;\mu\nu}(t_{i}+\tau,t_{i};\vec{p},\vec{q}_{1})L^{(f_{2})}_{P}(t_{i}-t_{P},\vec{p})\rangle\,,$
(19)

where translation invariance in the time direction has been used.

The computation of the pseudoscalar-vector two-point function (second diagram
in Fig. 2) is similar, except that we use a single stochastic source per time
slice, at the pseudoscalar insertion, but with support on the whole time
slice. We again use AMA to reduce the cost of the analysis. Since this
contribution is much smaller, we use only 16 sloppy inversions per exact solve
(2 for the strange quark). In this case, the final correlators for both
flavors $(f)$ in the PV loop are computed as

$C^{(3)}_{\mathrm{PV}_{f}-\mathrm{V};\mu\nu}(\tau,t_{P};\vec{p},\vec{q}_{1})=\frac{a}{T}\sum_{t_{i}}\
\langle
L^{(l-s)}_{V;\nu}(t_{i};\vec{q}_{1}-\vec{p})D^{(f)}_{V_{\mu}P}(t_{i}+\tau,t_{i}-t_{P};\vec{p},\vec{q}_{1})\\\
+L^{(l-s)}_{V;\mu}(t_{i}+\tau;-\vec{q}_{1})D^{(f)}_{V_{\nu}P}(t_{i},t_{i}-t_{P};\vec{p},\vec{p}-\vec{q}_{1})\rangle\,.$

The remaining correlation function (last diagram in Fig. 2) is the purely
disconnected contribution,
$C^{(3)}_{\mathrm{P}_{f}-\mathrm{V}-\mathrm{V};\mu\nu}(\tau,t_{P};\vec{p},\vec{q}_{1})$.
It can be easily computed from the loop functions in Eq. (18). Again, more
information is provided in Appendix B.

Finally, all correlation functions have been computed for two values of the
meson momentum ($\vec{p}=\vec{0}$ and $\vec{p}=(2\pi/L)\vec{z}$). The number
of photon momenta ($\vec{q}_{1}$) is tuned on each ensemble to efficiently
probe the ($Q_{1}^{2},Q_{2}^{2}$) plane up to 2 GeV2 for the pion and 4 GeV2
for the $\eta$ and $\eta^{\prime}$ mesons, see Fig. 3.

Figure 3: Kinematical reach for the pion, $\eta$ and $\eta^{\prime}$ meson at
the lattice spacing $a=0.0787\leavevmode\nobreak\ $fm. Each orbit corresponds
to a given value of $|\vec{q}_{1}|$ and $|\vec{q}_{2}|$. The whole orbit can
be generated by continuously varying the real parameter $\omega_{1}$. Black
and blue lines correspond to the pseudoscalar rest frame and moving frame
respectively.

##  IV Analysis

In Ref. [59], the amplitude $\widetilde{A}_{\mu\nu}(\tau)$ was written in
terms of a single scalar function $\widetilde{A}^{(1)}(\tau)$ and its
derivative

$\widetilde{A}_{\mu\nu}(\tau)=-iQ^{E}_{\mu\nu}\
\widetilde{A}^{(1)}(\tau)+P^{E}_{\mu\nu}\
\frac{\mathrm{d}\widetilde{A}^{(1)}}{\mathrm{d}\tau}(\tau)\,.$ (20)

The two coefficients $P^{E}_{\mu\nu}=iP_{\mu\nu}$ and
$Q^{E}_{\mu\nu}=(-i)^{n_{0}}Q_{\mu\nu}$ do not depend on $\omega_{1}$ and are
given by

$q_{\mu\nu}\equiv\epsilon_{\mu\nu\alpha\beta}q_{1}^{\alpha}q_{2}^{\beta}=P_{\mu\nu}\omega_{1}+Q_{\mu\nu}\,.$
(21)

In the analysis, all equivalent contributions to $\widetilde{A}^{(1)}$ are
averaged over and we explicitly use the Bose symmetry relation

$\widetilde{A}_{\mu\nu}(\tau;\vec{q}_{1},\vec{q}_{2})=\widetilde{A}_{\nu\mu}(-\tau;\vec{q}_{2},\vec{q}_{1})\,e^{-E_{P}\tau}\,.$
(22)

Once $\widetilde{A}^{(1)}(\tau)$, or equivalently
$\widetilde{A}_{\mu\nu}(\tau)$ is known, the integration over $\tau$ in Eq.
(3) is done using the Simpson rule. In Section IV.1, we will often display
individual contributions to $\widetilde{A}^{(1)}(\tau)$ that originate from
the different Wick contractions of Fig. 2. The scalar function
$\widetilde{A}^{(1)}(\tau)$ is strongly affected by staggered oscillations due
to the contribution of the parity partner state. We stress that these
oscillations do not disappear as we take the continuum limit but the
contribution of the parity partner state to the form factor does. Thus, for
clarity, we will often display the smeared function

$f(\tau)|_{\rm
smr}=\frac{1}{4}f(\tau-a)+\frac{1}{2}f(\tau)+\frac{1}{4}f(\tau+a)\,,$ (23)

instead of $f$ itself.

###  IV.1 Wick contractions

The contributions from the four Wick contractions of Fig. 2 to the amplitude
$\widetilde{A}^{(1)}(\tau)$ are shown in Fig. 4 for one ensemble at our finest
lattice spacing. In the isospin limit, the last two contractions, that involve
a pseudoscalar loop, vanish for the pion and contribute only to the $\eta$ and
$\eta^{\prime}$ transition form factors. In [59, 79], the authors have shown
that the second diagram, with a single vector loop, is numerically sub-
dominant. It is negative and contributes only at the percent level to the pion
TFF. In particular, this diagram would exactly vanish in the SU(3) flavor
limit. For the same reason, the last diagram of Fig. 2, that contains two
vector loops, is expected to give a very small contribution to the $\eta$ and
$\eta^{\prime}$ form factors. Thus, for the $\eta$ and $\eta^{\prime}$ TFFs,
most of the signal comes from the purely connected contribution and from the
disconnected contribution that contains a single pseudoscalar loop. Both
contributions have opposite sign and the cancellation is stronger for the
$\eta^{\prime}$.

Figure 4: Contributions from the four Wick contractions, displayed in Fig. 2,
to the integrand in Eq. (3), for the $\eta$ (left) and $\eta^{\prime}$ (right)
mesons. The figure corresponds to the choice $\omega_{1}=\omega_{s}\equiv
E_{P}/2$ that would lead to the double-virtual TFF with $Q_{1}^{2}=Q_{2}^{2}$.
The smearing procedure is defined by Eq. (23). The black circles represent the
sum over all contributions. Some contributions have been rescaled by a factor
50 for clarity. The results are given at our finest lattice spacing.

###  IV.2 Tail correction

Because of the finite time extent of the lattice and of the increase of the
statistical error, the integration over $|\tau|$, the time separation between
the two vector currents, is truncated at $\tau_{c}\approx 1.5$ fm. In Ref.
[59], the large-$\tau$ dependence of the amplitude was fitted assuming a
vector meson dominance (VMD) or a lowest meson dominance (LMD)
parametrization. This parametrization was then used to estimate the integrand
above the cut. With staggered quarks, this procedure is more complicated due
to the contribution of the parity partner state that appears with the usual
factor $(-1)^{\tau}$. This contribution, responsible for the presence of
oscillations at the level of the integrand, largely cancels between
consecutive time slices and vanishes in the TFFs in the continuum limit.
However, since we have little knowledge about the parity parter state
contribution, and to avoid complicated fits, we decided to follow a different
strategy. First, we compute the TFF neglecting the tail of the integrand above
the cut. The resulting form factor is fitted using an LMD parametrization in
the range $Q^{2}\leq Q^{2}_{\rm tail}=1$ GeV2 where the model provides an
acceptable description of our data. In a second step, this parametrization is
used to estimate the missing tail contribution at the level of the amplitude
(we use the equations derived in Appendix A of Ref. [59]). In principle, this
method can be applied iteratively until it converges. In practice, since we
use a conservative value of $\tau_{c}$, we find that at most two iterations
already suffice at our level of precision. To avoid any bias due to this tail
correction, we remove all points where the tail contributes more than 10% to
the integrand in Eq. (3).

###  IV.3 Finite-size effects

Figure 5: Pion TFF in the double-virtual regime using two different volumes
($L=3$ fm and $L=6$ fm), at our finest lattice spacing. Before (left) and
after (right) correcting for backward propagating pion.

The calculation of the disconnected diagram with a single pseudoscalar loop
(third diagram in Fig. 2) is numerically expensive and dominates the cost of
our simulations. In this section we study the possibility of using smaller
physical volumes to increase the statistical precision.

When using smaller volumes, two different kind of finite-size effects need to
be considered. First, it was observed in Ref. [58] that the finite-time extent
of the lattice potentially leads to large finite-time effects (FTE) due to
backward propagating pions. Fortunately, a simple procedure was proposed to
correct for this effect. Second, we need to consider finite-spatial-extent
effects. This can be done by looking at the pion TFF for which a high
statistical precision has been achieved on both large ($L=6$ fm) and small
($L=3$ fm) ensembles. Since the mass of the taste-singlet pseudoscalar meson
decreases as we approach the continuum limit, we focus on our finest lattice
spacing where the taste-singlet pion mass is approximately 200 MeV. The
results for $L=6$ fm and $L=3$ fm boxes are depicted on Fig. 5 before and
after FTE correction. In fact, when one corrects for FTE, the data for the two
different box sizes agree within statistical errors, even at the lowest
accessible virtualities $Q^{2}$. For the pion TFF, that is extracted
exclusively from $L\approx 6\leavevmode\nobreak\ $fm ensembles, finite-size
effects (FSE) are thus neglected.

Large volumes allow us to probe much smaller virtualities. This feature is
very important for the pion-pole contribution, but to a lesser extent for the
$\eta$\- and $\eta^{\prime}$-pole contributions. First
$a_{\mu}^{\mathrm{HLbL};\eta}$ and $a_{\mu}^{\mathrm{HLbL};\eta^{\prime}}$ are
less sensitive to the region close to the origin, see Table 2. Second, as the
pseudoscalar mass increases, lower virtualities can be reached, even in small
volumes. To estimate the FSE due to the use of small spatial volumes at fine
lattice spacings, the pion-pole contribution is computed from Eq. (16) but
using the weight functions of the $\eta^{(\prime)}$ meson. The relative
difference between small and large volumes is then taken as our systematic.

###  IV.4 Excited-state contributions

When using Eq. (7) to extract the $\eta$ form factor, instead of Eq. (11), the
$\eta^{\prime}$ is the first excited state and its relative contribution
decays exponentially with $\Delta
E\,t_{P}\equiv(E_{\eta^{\prime}}-E_{\eta})t_{P}$

$\widetilde{A}_{\mu\nu}^{(\eta);{\rm
eff}}(\tau,t_{P})=\widetilde{A}_{\mu\nu}^{(\eta)}(\tau)\left(1+\frac{E_{\eta}}{E_{\eta^{\prime}}}\frac{Z_{8}^{(\eta^{\prime})}}{Z_{8}^{(\eta)}}\frac{\widetilde{A}_{\mu\nu}^{(\eta^{\prime})}(\tau)}{\widetilde{A}_{\mu\nu}^{(\eta)}(\tau)}e^{-\Delta
E\,t_{P}}+\cdots\right)\,,$ (24)

where the ellipsis represents higher excited state contributions. The ratio of
overlap factors strongly depends on the choice of interpolating operator used
for the $\eta$ and, at the physical point, the mass gap is $\Delta E\approx
400$ MeV. With staggered quarks, because of taste breaking effects, the mass
gap between the two taste-singlet mesons decreases at large lattice spacings:
we find $\Delta E\approx 300$ MeV at our coarsest lattice spacing. With our
choice of interpolating operator, and at the physical $\eta^{(^{\prime})}$
masses, we find that $t_{P}\geq 1.0\leavevmode\nobreak\ $fm is needed to
reduce the exponentially suppressed term to the level of 10% and 1.35 fm for
5% (the target precision if one aims at 10% uncertainty on the pseudoscalar-
pole contribution to $a_{\mu}^{\rm hlbl}$). We conclude that large values of
$t_{P}$ are required unless the $\eta$ interpolating operator is carefully
chosen.

Despite the use of AMA and large statistics, the disconnected diagram that
contains a single pseudoscalar loop is noisy and becomes the main source of
statistical error. Thus, it is very useful to use as small values of $t_{P}$
as possible. If, instead of Eq. (7), one extracts the TFFs from a matrix of
correlation functions, the $\eta^{\prime}$ contribution to the $\eta$ TFF is
essentially suppressed (up to statistical precision in the determination of
the overlap factors). In Fig. 6, we compare the effective value of
$\widetilde{A}^{(1)}(\tau=0)$ as a function of $t_{P}$ using both methods. The
second approach, in black, where the $\eta^{\prime}$ contribution is
effectively removed, leads to a plateau at much earlier times. The value of
$t_{P}$, for each value of $\tau$, is chosen to minimize the statistical error
while keeping systematic errors under control (vertical line in Fig. 6). In
practice, $\widetilde{A}^{(1)}(\tau)$ it obtained from a smearing over three
values of $t_{P}$ using the same weights factors as in Eq. (23).

Figure 6: Effective value for the function $\widetilde{A}^{(1)}(\tau=0)$ as a
function of $t_{P}$ for the $\eta$ meson at the lattice spacing
$a=0.0787\leavevmode\nobreak\ $fm. The vertical green band corresponds to our
optimal choice. Black circles are obtained using the improved estimator given
by Eq. (11), while the blue squares are the corresponding data using the
estimator in Eq. (7). The orange line is the excited state prediction for the
estimator (7) as predicted by Eq. (8). Blue points are slightly shifted on the
right for clarity.

In Fig. 7, we present two stability plots where the pseudoscalar pole
contribution for the $\eta$ (black) and $\eta^{\prime}$ (red) are shown as a
function of $\Delta t_{P}$, a constant shift from our optimal values for
$t_{P}$ (negative values mean more agressive choices of $t_{P}$). With the
improved estimator, we hardly see any trend, except for the $\eta^{\prime}$ at
fine lattice spacing. For comparison, we also present the result for the
$\eta$ meson using the naive estimator (in gray). In this case, we observe a
clear slope due to excited-state contamination, and both methods to extract
$a_{\mu}^{\mathrm{HLbL};\eta}$ eventually agree at large $t_{P}$ although the
statistical error increases significantly.

Figure 7: Value of $a_{\mu}^{\mathrm{HLbL};P}$ as a function of $\Delta t_{P}$
using the improved estimator. The point $\Delta t_{P}=0$ corresponds to the
choice used in our analysis. Left: $a=0.1315$ fm. Right: $a=0.0787$ fm. The
open grey circles are the results for the $a_{\mu}^{\mathrm{HLbL};\eta}$ using
the estimator (7).

###  IV.5 Treatment of statistical and systematic errors

Statistical errors are estimated using the jackknife method after blocking.
All observables are computed on the same set of gauge field configurations and
using the same block size such that we can easily propagate statistical
errors. The lattice spacing was computed in [16], with a relative precision at
the permil level, but using a different blocking procedure. The lattice
spacing enters as a trivial factor in the conversion of the TFF and the
virtualities in physical units but also in the continuum extrapolations. We
have generated pseudo-samples with a gaussian distribution to propagate the
error on the lattice spacing. This error is always sub-dominant so that it is
safe to neglect the correlation of the lattice spacing with the quantities
studied here.

During the analysis several choices are made. First, to extract the form
factors on individual ensembles, we need to choose $\Delta t_{P}$ (for the
$\eta$ and $\eta^{\prime}$ mesons, see Section IV.4), the fitting range
$Q^{2}_{\rm tail}$ for the LMD or VMD model used for the tail correction, the
values of $\omega_{1}$ that are used to samples the orbits (see Eq. (3) and
above) and the strategy to compute the pseudoscalar energies and overlap
factors (see Section V.4.2). Second, in the continuum extrapolation, we can
choose the order of the $z$-expansion described in Section V.2, add higher
order terms in the continuum extrapolation or perform cuts in lattice
spacings. Several analyses have been performed by varying all these choices.
Since some fits do not include all correlations, the systematic error is
estimated assuming a flat weight among these representative variations. An
explicit example is given in Section V.4.2 where we compute the pseudoscalar-
pole contributions.

##  V Results

###  V.1 Transition form factors at a single lattice spacing

Results for the TFFs, at a single lattice spacing ($a=0.0787$ fm), are
presented in Fig. 8. For clarity, we present our results for two specific
kinematics: the single-virtual form factor with one real photon and the
double-virtual form factor where both photons share the same virtuality
$Q_{1}^{2}=Q_{2}^{2}$. In both cases, we observe a good agreement between the
pseudoscalar rest frame and moving frame. In [59] the motivation for adding a
non-vanishing pseudoscalar momentum was to better probe the
($Q_{1}^{2},Q_{2}^{2}$) plane. This is crucial for the pion because of its low
mass (left panel in Fig. 3). For the heavier $\eta$ and $\eta^{\prime}$, a
good coverage of the ($Q_{1}^{2},Q_{2}^{2}$) plane, even for the single-
virtual TFF, is already reached in the meson rest frame.

For the pion, the signal essentially comes from the fully connected
contribution with a small negative contribution, at the percent level, from
the second disconnected diagram in Fig. 2 (see Appendix B). Thus, a very high
statistical precision can be achieved. Nevertheless, the signal-to-noise ratio
deteriorates at low virtualities. This is of particular relevance when
computing the pseudoscalar-pole contribution to $a^{\rm HLbL}_{\mu}$ for which
the low-Q2 region dominates. This point is further discussed in Section V.4.
The signal for the $\eta$ and $\eta^{\prime}$ TFFs is noisier for two reasons.
First, the disconnected diagram with a single pseudoscalar loop (third diagram
in Fig. 2) is large and enters with an opposite sign as compared to the fully
connected contribution. Second, the overlap factors and the energies of the
pseudoscalar mesons also contribute to the error.

In the following section, we extrapolate our data to the continuum limit.

Figure 8: TFF of the $\pi^{0}$ (top panel), $\eta$ (second panel) and
$\eta^{\prime}$ (bottom panel) in the double (left) and single (right) virtual
regimes. Black filled circles and blue filled circles indicate respectively
$\vec{p}=\vec{0}$ and $\vec{p}=\frac{2\pi}{L}(0,0,1)$. The blue bands are fit
based on the parametrization given by Eq. (25) at a single lattice spacing
($a=0.0787\leavevmode\nobreak\ $fm).

###  V.2 Parametrization of ${\cal F}_{P\gamma^{*}\gamma^{*}}$ and
extrapolation to the physical point

To extrapolate the form factor to the continuum limit, we fit our lattice data
using the modified $z$-expansion introduced in [59, 80]

$P(Q_{1}^{2},Q_{2}^{2})\ {\cal
F}_{P\gamma^{*}\gamma^{*}}(-Q_{1}^{2},-Q_{2}^{2})=\\\
\sum_{n,m=0}^{N}c_{nm}(a)\,\left(z_{1}^{n}+(-1)^{N+n}\frac{n}{N+1}\,z_{1}^{N+1}\right)\,\left(z_{2}^{m}+(-1)^{N+m}\frac{m}{N+1}\,z_{2}^{N+1}\right)\,,$
(25)

where $z_{k}$ are conformal variables

$\displaystyle
z_{k}=\frac{\sqrt{t_{c}+Q_{k}^{2}}-\sqrt{t_{c}-t_{0}\vphantom{Q_{k}^{2}}}}{\sqrt{t_{c}+Q_{k}^{2}}+\sqrt{t_{c}-t_{0}\vphantom{Q_{k}^{2}}}},\quad
k=1,2,$ (26)

$c_{nm}$ are coefficients symmetric in $m$ and $n$, $t_{c}=4m_{\pi}^{2}$ maps
the branch cut of the TFF onto the unit circle $\left|z_{k}\right|=1$ and
$t_{0}$ is a free parameter. The choice
$t_{0}=t_{c}(1-\sqrt{1+Q_{\max}^{2}/t_{c}})$ ensures that the maximum value of
$\left|z_{k}\right|$ is minimized in the momentum range $Q^{2}<Q_{\max}^{2}$.
Finally, the function $P(Q_{1}^{2},Q_{2}^{2})$ is an arbitrary analytic
function. In [59] it was argued that a convenient choice is

$P(Q_{1}^{2},Q_{2}^{2})=1+\frac{Q_{1}^{2}+Q_{2}^{2}}{\Lambda^{2}}\,,$ (27)

such that the form factor parametrization decreases as $1/Q^{2}$ at large
virtualities, even at finite value of $N$, as expected from the operator
product expansion (OPE) [81, 82, 83, 84, 85]. For the pion and the $\eta$
TFFs, $\Lambda=775\leavevmode\nobreak\ $MeV is set to the $\rho$ meson mass.
For the $\eta^{\prime}$, we use $\Lambda=1\leavevmode\nobreak\ $GeV. In
practice, we do not observe a strong dependence on this parameter.

With staggered quarks, we expect discretization errors to enter quadratically
[70, 71]. They are taken into account by expanding the $c_{nm}$ coefficients
as

$c_{nm}(a)=c_{nm}(0)\left(1+\gamma_{nm}(a\Lambda_{\rm QCD})^{2}\right)\,.$
(28)

with $\Lambda_{\rm QCD}=0.5\leavevmode\nobreak\ \mathrm{GeV}$ taken as a
typical QCD scale. The data for the pion is the most precise. In this case, to
check for possible higher order discretization effects and to estimate the
associated systematic error, additional fits that include a quartic term
$\delta_{nm}a^{4}$ have been performed. However, to avoid over-fitting, those
additional parameters are included only when the error on the corresponding
quadratic coefficient, $\gamma_{nm}$, is well below 50%. For the pion, the
data is also limited to virtualities $Q^{2}\lesssim 2.2\leavevmode\nobreak\
\mathrm{GeV}^{2}$. To stabilize the fit at large virtualities, the additional
constraint $Q^{2}{\cal
F}_{\pi^{0}\gamma^{*}\gamma^{*}}(-Q^{2},-Q^{2})=0.062(15)\leavevmode\nobreak\
$GeV at $Q^{2}=20\leavevmode\nobreak\ \mathrm{GeV}^{2}$ is imposed. We believe
that this 25% uncertainty is conservative as it is larger than the difference
between the TFF value at $Q_{1}^{2}=Q_{2}^{2}=2\leavevmode\nobreak\
\mathrm{GeV}^{2}$ and the asymptotic value [81, 82, 83]

$Q^{2}{\cal
F}_{\pi^{0}\gamma^{*}\gamma^{*}}(-Q^{2},-Q^{2})\xrightarrow[Q^{2}\to\infty]{}\frac{2F_{\pi}}{3}$
(29)

predicted by the OPE. For the pion decay constant in the chiral limit, we use
$F_{\pi}=92.4\leavevmode\nobreak\ $MeV.

The main analysis is done using $N=2,3$ for the pion and for the $\eta$ meson,
and with $N=1,2$ for the less precise $\eta^{\prime}$ TFF. As can be seen in
Fig. 8, the number of points included in the fits is very large. Furthermore,
many ensembles are fitted simultaneously. It is clear that a fully correlated
fit is not feasible. Thus, we decided to perform uncorrelated fits but to
estimate the error through the jackknife procedure. This strategy prevents us
from giving a clear meaning to the fit quality, but leads to reliable error
estimates [86, 87].

The final results are depicted in Fig. 9. The fit parameters corresponding to
one of our systematic variations (see Section IV.5), and the associated
correlation matrix, for each pseudoscalar, are given in Appendix A. The
uncorrelated chi-squared per degree of freedom of the global fits are
typically 1.15, 0.95 and 0.85 for the $\pi^{0}$, $\eta$ and $\eta^{\prime}$
cases. We note that $\chi^{2}/\mathrm{d.o.f.}\approx 1$ are still expected as
the global fits include uncorrelated ensembles. In Fig. 10, we show the
continuum extrapolation of the TFFs for three different sets of virtualities
$(Q_{1}^{2},Q_{2}^{2})$. In this case, each ensemble is fitted individually
using Eq. (25) and the extrapolation is compared with the global fit
procedure. In all cases, we observe a quadratic scaling in the lattice
spacing. The normalization of the TFFs are discussed in the following section.

Figure 9: TFFs at the physical point using the $z-$expansion from Eq. (25).
Errors are statistical only. The horizontal black line corresponds to the OPE
prediction. The dispersive result is extracted from [61], the Canterbury
appoximant (CA) result from [88] and the Dyson-Schwinger equation (DSE) result
comes from [89]. Measurements from CELLO [50], CLEO [51] and L3 [54] are shown
for comparison in the single-virtual case. The preliminary results from the
BES III experiment have been extracted from [3]

In the single-virtual case, our result for the $\eta$ TFF shows a slight
tension with the smallest CELLO bin. This feature is reflected by the tension
in the normalization for the TFF discussed in the next section. For the
$\eta^{\prime}$ meson, our result is slightly above the CELLO data although
the tension is well below two combined standard deviations. The pion TFF is in
very good agreement with experimental measurements. For the pion and
$\eta^{\prime}$ TFFs we find reasonable agreement with previous determinations
based on lattice calculations of the pion TFF [59], the dispersive framework
[61], Canterbury approximants [88] and Dyson-Schwinger equations [89]. In the
case of the $\eta$ TFF, we observe a tension below $0.4\leavevmode\nobreak\
\mathrm{GeV}^{2}$ that is explained by our small value of two-photon decay
rate as compared to experiments.

Figure 10: Continuum extrapolations of the TFFs for three kinematics
$(Q_{1}^{2},Q_{2}^{2})$. Black and red points are obtained by fitting
individual ensembles using Eq. (25). The blue points are the result of the
global fit procedure described in the text. Errors are statistical only.

###  V.3 Normalization of the form factors

When both photons are real, the form factor is related to the decay rate of
the pseudoscalar meson to two photons. At leading order in QED

$\Gamma(P\to\gamma\gamma)=\frac{\pi\alpha^{2}m_{P}^{3}}{4}\mathcal{F}^{2}_{P\gamma\gamma}(0,0)$
(30)

where $\alpha$ is the fine-structure constant. The neutral pion decay width
has been measured by the PrimEx [90] and PrimEx-II [57] experiments with a
precision of 1.5% leading to the combined result
$\Gamma(\pi^{0}\to\gamma\gamma)=7.802(52)_{\rm stat}(105)_{\rm syst}$ eV. It
translates to a relative precision of 0.75% on the normalization of the form
factor. From the PDG average [91], the $\eta$ and $\eta^{\prime}$ decay widths
are $\Gamma(\eta\to\gamma\gamma)=516\pm 18$ eV and
$\Gamma(\eta^{\prime}\to\gamma\gamma)=4.28\pm 0.19$ keV. It corresponds to a
relative precision of 1.7% and 2.2% respectively. In the case of the $\eta$
meson, a significant tension exists between the Primakoff measurement [92] and
the results from collider experiments [93, 94, 95, 96, 97]. A new measurement
by PrimEx-Eta is expected in the comings years [98].

Figure 11: Continuum extrapolation of the normalization of the form factor
${\cal F}_{P\gamma^{*}\gamma^{*}}(0,0)$ from one typical variation described
in Section IV.5. The blue points are the results obtained from the global fit
strategy. Black and red points are obtained by fitting the form factor on each
ensemble separately. The corresponding extrapolation is then shown in green.
For the $\eta$ and $\eta^{\prime}$ mesons, the red points correspond to large-
volume ensembles only.

From our global fit procedure, described in Section V.2, we obtain

$\displaystyle\Gamma(\pi^{0}\to\gamma\gamma)$
$\displaystyle=7.11\,\pm\,0.44_{\rm stat}\,\pm\,0.21_{\rm
syst}\leavevmode\nobreak\ \mathrm{eV}$ (31a) $\displaystyle\Gamma(\eta^{\
}\to\gamma\gamma)$ $\displaystyle=338\,\pm\,94_{\rm stat}\,\pm\,35_{\rm
syst}\leavevmode\nobreak\ \mathrm{eV}$ (31b)
$\displaystyle\Gamma(\eta^{\prime}\to\gamma\gamma)$
$\displaystyle=3.4\,\pm\,1.0_{\rm stat}\,\pm\,0.4_{\rm
syst}\leavevmode\nobreak\ \mathrm{keV}\,.$ (31c)

Our results for the pion and the $\eta^{\prime}$ mesons are compatible with
the experimental values at the level of 1.4 and 0.8 standard deviation
respectively. For the $\eta$ meson, our value presents a slight tension at the
level of 1.7 standard deviations and is smaller in magnitude. The impact on
the pseudoscalar pole-contribution will be discussed in the following section.

We note that, at our finer lattice spacings, our results are based on
relatively small physical volumes. Furthermore, FSE are expected to be the
largest at small virtualities and might affect the normalization of the form
factor more strongly. In addition to the global fit, we have also performed
individual fits on each ensemble separately. The large volume ensembles, shown
in red in Fig. 11, fully support the small decay rate for the $\eta$ meson.
Furthermore, we remind the reader that we do not see any significant FSE for
the pion where the calculation has been performed on both large and small
volumes at all lattice spacings. We also tried to fit our data using the LMD
model described in [59] and we obtain a fully compatible result,
$\Gamma(\eta^{\ }\to\gamma\gamma)=301\,\pm\,50_{\rm stat}$ eV. Finally, we
note that the ETM result $\Gamma(\eta^{\ }\to\gamma\gamma)=323\pm
88\leavevmode\nobreak\ \mathrm{eV}$ [62] is also smaller than the experimental
value based on collider experiments. However, this result is based on a single
ensemble and it is difficult to predict the size of the discretization
effects.

###  V.4 Pseudoscalar-pole contribution to HLbL scattering in $(g-2)_{\mu}$

The master formula Eq. (16) requires the knowledge of the TFF at any spacelike
virtuality and the dominant contribution comes from the low-energy region with
$Q^{2}\lesssim 2\leavevmode\nobreak\ \mathrm{GeV}^{2}$. Since the single-
virtual TFF appears explicitly, this specific kinematic is important. For the
pion, it was noticed in [59] that the single-virtual TFF is more difficult to
extract on the lattice and this observation motivated the inclusion of a
moving frame for the pion to improve the covering of large virtualities. The
situation tends to be better for the $\eta$ and $\eta^{\prime}$ mesons due to
their heavier masses: we probe larger virtualities and the signal deteriorates
more slowly at large virtualities.

A detailed analysis of the weight functions $w_{1}$ and $w_{2}$ has been done
in [99] and we simply mention their most relevant properties. First, the
second weight function $w_{2}$ is typically one order of magnitude smaller
than $w_{1}$ and, at large virtualities, the weight function $w_{1}$ goes to
zero as

$w_{1}(Q_{1},Q_{2},\tau)\underset{Q_{1}^{2}\to\infty}{\sim}Q_{1}^{-1}\,,\quad
w_{1}(Q_{1},Q_{2},\tau)\underset{Q_{2}^{2}\to\infty}{\sim}Q_{2}^{-2}\,,\quad
w_{1}(Q,Q,\tau)\underset{Q^{2}\to\infty}{\sim}Q^{-2}\,.$ (32)

In particular, there is a long tail at constant value of $Q_{2}^{2}$. Second,
the only dependence on the pseudoscalar mass appears as pre-factors
$(Q_{2}^{2}+m_{P}^{2})^{-1}$ and $((Q_{1}+Q_{2})^{2}+m_{P}^{2})^{-1}$ in
$w_{1}$ and $w_{2}$ respectively. Thus, we expect the saturation of the
integrand, as a function of the $Q_{1}$ and $Q_{2}$, to be slower for heavier
pseudoscalar mesons. In [99], it was noted that, for the bulk of the
contribution, the relation $w_{1}^{\eta}=(1/6)w_{1}^{\pi}$ and
$w_{1}^{\eta^{\prime}}=(1/2.5)w_{1}^{\eta}$ approximately holds.

####  V.4.1 Lattice result

Figure 12: Continuum extrapolation of the pseudoscalar-pole contribution
$a_{\mu}^{\mathrm{HLbL};P}$ from one typical analysis variation described in
Section IV.5. The black and red points are obtained by fitting the TFFs for a
single ensemble to the parametrization of Eq. (25). For the $\eta$ and
$\eta^{\prime}$ mesons, the red points correspond to large-volume ensembles
only. We emphasize that single-ensemble fits tend to be instable at large
$Q^{2}$, especially for the pion. Thus the global fit strategy is always our
preferred choice. It leads to the blue points at $a=0$

We follow two strategies to perform the continuum extrapolation. In the first
strategy, we start with computing $a_{\mu}^{\mathrm{HLbL};P}$ ensemble per
ensemble and then extrapolate the observable to the continuum limit. The
result is displayed by the green band in Fig. 12. The alternative strategy is
to evaluate Eq. (16) using the parametrization of the TFF in the continuum
limit. The result is shown in blue in Fig. 12. Both methods lead to compatible
results and we quote the second method as our preferred strategy because the
fits involved are more stable. It leads to

$\displaystyle a_{\mu}^{\mathrm{HLbL};\pi^{0}}$
$\displaystyle=(57.8\,\pm\,1.8_{\rm stat}\,\pm\,0.9_{\rm syst})\times
10^{-11},$ (33a) $\displaystyle a_{\mu}^{\mathrm{HLbL};\eta}\,$
$\displaystyle=(11.6\,\pm\,1.6_{\rm stat}\,\pm\,0.5_{\rm syst}\,\pm\,1.1_{\rm
FSE})\times 10^{-11},$ (33b) $\displaystyle
a_{\mu}^{\mathrm{HLbL};\eta^{\prime}}$ $\displaystyle=(15.7\,\pm\,3.9_{\rm
stat}\,\pm\,1.1_{\rm syst}\,\pm\,1.3_{\rm FSE})\times 10^{-11}.$ (33c)

Summing all contributions, and taking correlation into account, we obtain the
pseudoscalar-pole contribution

$a_{\mu}^{\mathrm{HLbL};P}=(85.1\,\pm 4.7_{\rm stat}\,\pm\,2.3)\times
10^{-11}\,.$ (34)

We observe small anti-correlations (with a correlation coefficient of order
$-0.1$) between the $\pi^{0}$ and $\eta$-pole contributions as well as between
the $\pi^{0}$ and $\eta^{\prime}$-pole contributions. The correlation
coefficient between the $\eta$ and $\eta^{\prime}$-pole contribution is 0.35.
The pion-pole contribution is in good agreement with previous lattice
calculations [58, 59] that found
$a_{\mu}^{\mathrm{HLbL};\pi^{0}}=59.7(3.6)\times 10^{-11}$. It is compatible
at the level of 1.5$\sigma$ with the data-driven estimate [61],
$a_{\mu}^{\mathrm{HLbL};\pi^{0}}=(62.6^{+3.0}_{-2.5})\times 10^{-11}$ used in
the average presented in [3]. The $\eta$-pole contribution is in slight
tension with the White Paper estimate based on Canterbury approximants [88],
$a_{\mu}^{\mathrm{HLbL};\eta}=16.3(1.4)\times 10^{-11}$, but also with the
results based of Dyson-Schwinger equations:
$a_{\mu}^{\mathrm{HLbL};\eta}=15.8(1.2)\times 10^{-11}$ [89] and
$a_{\mu}^{\mathrm{HLbL};\eta}=14.7(1.9)\times 10^{-11}$ [100]. This tension is
mostly due to the momentum region below 0.5 $\mathrm{GeV}$ as indicated by the
low value of the decay rate. Finally, for the $\eta^{\prime}$-pole, our result
is in good agreement with the Canterbury appoximant (CA) estimate [88]:
$a_{\mu}^{\mathrm{HLbL};\eta^{\prime}}=14.5(1.9)\times 10^{-11}$ and with the
result based on Dyson-Schwinger equations:
$a_{\mu}^{\mathrm{HLbL};\eta^{\prime}}=13.3(0.9)\times 10^{-11}$ [89] and
$a_{\mu}^{\mathrm{HLbL};\eta^{\prime}}=13.6(0.8)\times 10^{-11}$ [100]. A
summary plot is shown in Fig. 14.

In Table 2 we provide the relative contribution of the pseudoscalar-pole
contributions $a_{\mu}^{\rm hlbl,P}$ as a function of the momentum cutoff
$Q_{\rm cut}$. Here, the TFF is assumed to be zero if either $|Q_{1}|>Q_{\rm
cut}$ or $|Q_{2}|>Q_{\rm cut}$. Note that it differs from a cut in the
integration range in the master equation (16) as one of the arguments of one
of the TFFs in each term of the sum is $Q_{1}^{2}+Q_{2}^{2}$. Since the TFF is
the central quantity computed in this work, we believe that this choice is
more meaningful. If the contribution from the region
$|Q|<0.25\leavevmode\nobreak\ $GeV represents more than 50% for the pion, it
drops to 31% for the $\eta$ and even 20% for the $\eta^{\prime}$. Concerning
the contribution to the total error, we observe a similar scaling with $Q_{\rm
cut}$. In particular, for the pion, half of the error originates from the
region $|Q|<0.25\leavevmode\nobreak\ $GeV where the signal to noise ratio
deteriorates.

Table 2: Relative contribution of the pseudoscalar-pole contributions $a_{\mu}^{\rm hlbl,P}$ as a function of the momentum cutoff $Q_{\rm cut}$ (see the main text for explanations).  $Q_{\rm cut}$ [GeV] |  $\pi^{0}$ |  $\eta$ | $\eta^{\prime}$
---|---|---|---
0.25 | 54% | 31% | 20%
0.50 | 71% | 50% | 38%
0.75 | 80% | 61% | 51%
1.00 | 84% | 68% | 60%
1.50 | 89% | 76% | 71%
2.00 | 92% | 81% | 78%
3.00 | 94% | 86% | 85%
5.00 | 97% | 91% | 92%

####  V.4.2 Discussion of systematic errors

In this section, we discuss the different analysis variations that have been
included to estimate the systematic errors using the procedure presented in
Section IV.5. A summary plot is provided in Fig. 13.

For the pion, three different procedures have been used to extract the mass
and overlap factor with our interpolating operator: we start with a two-
exponential fit in a time range that leads to a good chi-square (2exp). In a
second step, this fit is used to determine a fit range where the excited state
contribution is negligible and the data are fitted by a single-exponential
(1exp). Finally, we perform a combined fit where both pion momenta are fitted
simultaneously assuming the validity of the dispersion relation (comb).
Similarly, for the $\eta^{(\prime)}$ mesons, three different procedures,
described in [63], have been considered. In this case, the masses and overlap
factors are obtained by fitting a correlation matrix with both vanishing and
non-vanishing pseudoscalar momenta. At vanishing momentum, it is often
advantageous to consider the forward derivative $D(t)=C(t+a)-C(t)$ that
significantly reduces the statistical noise. Thus, the three variations refer
to the combinations $(C,C)$, $(D,C)$ and $(D,D)$ where $(D,C)$ means that the
forward derivative $D$ is used at vanishing momenta while the original
correlator $C$ is used for the moving frame.

For all TFFs, three samplings of the orbit shown in Fig. 3 have been
considered. They differ by the number of $\omega$-values for each orbit,
$N_{\omega}=15,20,25$. We also considered three methods for the tail
correction, that differ by the value of $Q^{2}_{\rm
tail}=(0.75,0.75,1.25)\leavevmode\nobreak\ $GeV and/or the model (VMD, LMD,
LMD) used to fit the lattice data in the region
$Q_{1}^{2},Q_{2}^{2}<Q^{2}_{\rm tail}$.

For the $\eta$ and $\eta^{\prime}$ TFFs, we also considered two different
values of $\Delta t_{P}$. Our optimal choice, from Section IV.4, and a more
conservative choice.

In the continuum extrapolation of the pion TFFs, we have considered $a^{4}$
correction as explained in Section V.2. For the $\eta$ and $\eta^{\prime}$
mesons, the lattice data are not precise enough to test higher order terms in
the continuum extrapolation. Instead, we have performed cuts in the lattice
spacing by removing all points with $a>0.1\leavevmode\nobreak\ $fm.

The impact of each particular choice is shown in Fig. 13 where each point is
obtained by averaging over all variations that are compatible with this
choice. Our total error is dominated by statistics, especially for the $\eta$
and $\eta^{\prime}$ mesons.

Figure 13: Study of the systematic error on $a_{\mu}^{\mathrm{HLbL};P}$ for
the pion, $\eta$ and $\eta^{\prime}$. The orange, blue and grey bands
correspond to the systematic, statistical and total error respectively. For
each point, the result is obtained from a model average over all variations
that contain the specific choice.

##  VI Conclusion

Figure 14: Overview of recent estimates for the pion, $\eta$ and
$\eta^{\prime}$-pole contribution to HLbL scattering in muon $g-2$ in units of
$10^{-11}$. The data are extracted from the Mainz lattice result [59], the
data-driven approach (Disp 18) [61], the Canterbury approximants estimate (CA
19) [88], Dyson Schwinger estimates (DS 19) [89] and (DS 20) [100]. The last
line is the 2020 White Paper estimate [3].

In this paper, we have presented a first ab-initio calculation of the $\eta$
and $\eta^{\prime}$ transition form factors at the physical point, in the
kinematical range relevant for computing the HLbL contribution to the
anomalous magnetic moment of the muon. Results for the pion TFF are also
provided. As compared to the pion TFF, the extraction of the $\eta$ and
$\eta^{\prime}$ TFFs requires a proper treatment of the mixing in the
$\eta-\eta^{\prime}$ system [63]. Thanks to the implementation of advanced
numerical methods, we have been able to compute and get a signal for all
disconnected contributions: the diagram with a single pseudoscalar loop is by
far dominant, the other contributions are SU(3)f suppressed and contribute at
the percent level.

Our result for the pion TFF is in good agreement with previous lattice
calculations [58, 59] and with experimental data in the single-virtual regime.
Concerning the $\eta$ TFF, our result turns out to be in good agreement with
previous estimates based on Canterbury approximants [88] or Dyson-Schwinger
equations [89, 100] for virtualities above $0.4\leavevmode\nobreak\
\mathrm{GeV}^{2}$. However, we find a slight tension, at the level of
1.7$\sigma$, for the normalization of the $\eta$ TFF. In the current setup,
extracting the normalization of the TFF is challenging as it requires an
extrapolation with data becoming noisier as we approach vanishing
virtualities. Additional larger-volume ensembles could further improve the
result. The single-virtual $\eta$ form factor is also compatible with
experimental measurements by CELLO [50] and CLEO [51], with a slight tension
for the lowest CELLO bin. Finally, the result for the $\eta^{\prime}$ TFF is
in good agreement with other approaches in the double-vritual regime. In the
single-virtual regime, our result is compatible (at the level of 1-2$\sigma$)
with existing measurements from CELLO [50], CLEO [51], BaBar [55, 53] and the
L3 experiment at LEP [54]. We observe a slight tension with the results based
on the Dyson-Schwinger equation [89].

Using our parametrization of the TFFs at the physical point, we provide an ab-
initio estimate of the light-pseudoscalar meson-pole contribution to HLbL. For
the pion, we confirm the previous lattice calculations with Wilson-clover
quarks [59] where a good agreement was observed with the data-driven estimate
of Ref. [61]. The extraction of the $\eta$ and $\eta^{\prime}$ pole
contributions is much more challenging but a precision of 17% and 27% has been
achieved. Previous estimates of the $\eta^{(^{\prime})}$ pseudoscalar-pole
were based on Dyson-Schwinger equations [89, 100] or Canterbury approximants
[88], relying on the availability of experimental data. For the $\eta$ meson
we find a slightly smaller value, although compatible at the level of 1.9
combined standard deviations. This can be explained by our lower value for the
normalization of the form factor as compared to the experimental value quoted
in [56].

Summing all contributions in Eq. (33), we obtain the main result of this paper

$a_{\mu}^{\mathrm{HLbL};P}=(85.1\,\pm 4.7_{\rm stat}\,\pm\,2.3)\times
10^{-11}\,.$ (35)

This value is compatible with the White-Paper estimate,
$a_{\mu}^{\mathrm{HLbL};P}=93.8^{+4.0}_{-3.6}\times 10^{-11}$ [3] at the level
of 1.4 combined standard deviations, and with a competitive uncertainty.

In the future, adding large-volume ensembles at smaller lattice spacings for
the $\eta$ and $\eta^{\prime}$ TFFs would help to better constrain the
normalization of the form factor and to clarify the tension on the
$\eta\to\gamma\gamma$ decay rate. The main advantage of large physical volume
is to provide a better sampling of the low-virtuality region close to the
origin. As shown on Table 2, reducing the error on the decay rate has a strong
impact on the $\pi^{0}$ and $\eta$-pole contributions and, to a lesser extent,
on the $\eta^{\prime}$-pole contribution. Another possibility would be to use
the framework presented in [101, 102] to directly compute the normalization of
the TFF.

###### Acknowledgements.

We thank all the members of the Budapest-Marseile-Wuppertal collaboration for
helpful discussions and the access to the gauge ensembles used in this work.
This publication received funding from the Excellence Initiative of Aix-
Marseille University - A*Midex, a French “Investissements d’Avenir” programme,
AMX-18-ACE-005 and from the French National Research Agency under the contract
ANR-20-CE31-0016. The computations were performed on Joliot-Curie at CEA’s
TGCC, on Jean Zay at IDRIS, on SuperMUC-NG at Leibniz Supercomputing Centre in
München, on HAWK at the High Performance Computing Center in Stuttgart and on
JUWELS at Forschungszentrum Jülich. We thank GENCI (grants A0080511504,
A0100511504 and A0120511504) and the Gauss Centre for Supercomputing (projects
pn73xi and wprecision) for awarding us computer time on these machines. Centre
de Calcul Intensif d’Aix-Marseille (CCIAM) is acknowledged for granting access
to its high performance computing resources.

##  Appendix A Coefficients of the $z$-expansion

In this appendix we provide the coefficients of the $z$-expansion and the
associated correlation matrix as described in Section V.2. The results
correspond to one of the many analysis variations described in Section IV.5
and the errors are statistical only.

Table 3: Coefficients of the $z$-expansion at the physical point for the pion, $\eta$ and $\eta^{\prime}$ mesons with $N=2$. We also quote the value of $Q_{\max}^{2}$ and $\Lambda$ that enters in Eq. (25). Meson | $c_{00}$ | $c_{01}$ | $c_{11}$ | $c_{20}$ | $c_{21}$ | $c_{22}$ | $Q_{\max}^{2}$ | $\Lambda$
---|---|---|---|---|---|---|---|---
$\pi^{0}$ | $+0.2372(22)$ | $-0.064(05)$ | $-0.26(11)$ | $+0.082(55)$ | $+0.040(65)$ | $-0.060(53)$ | 2.2 GeV2 | 775 MeV
$\eta$ | $+0.198(11)$ | $-0.021(15)$ | $-0.39(15)$ | $+0.154(94)$ | $-0.22(19)$ | $+0.37(47)$ | 5 GeV2 | 775 MeV
$\eta^{\prime}$ | $+0.296(31)$ | $-0.055(25)$ | $-0.29(16)$ | $+0.001(103)$ | $+0.04(26)$ | $-0.01(71)$ | 5 GeV2 | 1 GeV

$\mathrm{cor}(c_{nm})_{\pi^{0}}=\begin{pmatrix}+1.00&-0.51&+0.20&-0.26&+0.37&-0.27\\\
-0.51&+1.00&-0.42&+0.40&-0.61&+0.45\\\ +0.20&-0.42&+1.00&-0.97&+0.70&+0.20\\\
-0.26&+0.40&-0.97&+1.00&-0.81&-0.10\\\ +0.37&-0.61&+0.70&-0.81&+1.00&-0.43\\\
-0.27&+0.45&+0.20&-0.10&-0.43&+1.00\\\ \end{pmatrix}\ $ (36)

$\mathrm{cor}(c_{nm})_{\eta}=\begin{pmatrix}+1.00&+0.08&+0.37&-0.42&-0.04&+0.19\\\
+0.08&+1.00&-0.25&+0.23&-0.74&+0.36\\\ +0.37&-0.25&+1.00&-0.89&+0.40&-0.02\\\
-0.42&+0.23&-0.89&+1.00&-0.44&-0.21\\\ -0.04&-0.74&+0.40&-0.44&+1.00&-0.60\\\
+0.19&+0.36&-0.02&-0.21&-0.60&+1.00\\\ \end{pmatrix}$ (37)

$\mathrm{cor}(c_{nm})_{\eta^{\prime}}=\begin{pmatrix}+1.00&-0.06&+0.01&-0.20&-0.33&+0.21\\\
-0.06&+1.00&-0.51&+0.36&-0.63&+0.17\\\ +0.01&-0.51&+1.00&-0.78&+0.64&-0.22\\\
-0.20&+0.36&-0.78&+1.00&-0.38&-0.29\\\ -0.33&-0.63&+0.64&-0.38&+1.00&-0.54\\\
+0.21&+0.17&-0.22&-0.29&-0.54&+1.00\\\ \end{pmatrix}$ (38)

##  Appendix B Disconnected contributions : estimator for the vector loops

In this appendix we compare different estimators to evaluate the vector loops
that appear in the second and fourth disconnected diagrams of Fig. 2. We also
provide estimates of those sub-leading diagrams to the pion-pole contribution
to $a_{\mu}^{\rm hlbl}$.

###  B.1 Noise reduction techniques

In the isospin limit, we are interested in the vector loop function for the
light-minus-strange contribution

$\displaystyle L^{(l-s)}_{V;\mu}(x_{0};\vec{q})$
$\displaystyle=\sum_{\vec{x}}\mathrm{Tr}\left[\gamma_{\mu}D^{-1}_{l}(x,x)-\gamma_{\mu}D^{-1}_{s}(x,x)\right]e^{i\vec{q}\cdot\vec{x}}$
(39)
$\displaystyle=(m_{s}-m_{l})\sum_{\vec{x}}\mathrm{Tr}\left[\gamma_{\mu}D^{-1}_{l}D^{-1}_{s}(x,x)\right]e^{i\vec{q}\cdot\vec{x}}$
(40)

where $D_{f}$ denotes the Dirac operator for a fermion of mass $m_{f}$. In
this appendix we use Wilson-like fermion notations. A standard estimator is to
use the same set of $N$ volume stochastic sources for both light and strange
contributions [103]

$L^{(l-s)}_{V;\mu}(x_{0};\vec{q})=(m_{s}-m_{l})\frac{1}{N}\sum_{n=1}^{N}\sum_{\vec{x}}\left[\eta_{n}^{{\dagger}}(x)\gamma_{\mu}D^{-1}_{l}D^{-1}_{s}\eta_{n}(x)\right]e^{i\vec{q}\cdot\vec{x}}$
(41)

where the volume sources satisfy
$\langle\eta(x)\eta^{{\dagger}}(y)\rangle=\delta_{x,y}$. This estimator is
noisy and typically thousands of inversions are required to reach the gauge
noise. In [78], the authors introduced the split-even stochastic estimator

$L^{(l-s)}_{V;\mu}(x_{0};\vec{q})=(m_{s}-m_{l})\frac{1}{N}\sum_{n=1}^{N}\sum_{\vec{x}}\left[\left(\eta_{n}^{{\dagger}}D_{l}^{-1}\right)(x)\gamma_{\mu}\left(D^{-1}_{s}\eta_{n}\right)(x)\right]e^{i\vec{q}\cdot\vec{x}}$
(42)

which drastically reduces the variance of the estimator for a comparable
numerical cost. Finally, both estimators can be combined with LMA where the
low-mode contribution is computed exactly, without adding stochastic noise.
When the low modes of the Dirac operator are used to deflate the lattice Dirac
operator, the overhead of LMA is small.

In this appendix, we compare the performances of the four estimators discussed
above. We use the variance of the loop function as a benchmark of the
estimators. Since the error on the loop function is constant in time, we
average the variance over all time slices and all directions $\mu=0,1,2,3$.
The results are obtained on the $L=6\leavevmode\nobreak\ $fm ensemble at a
lattice spacing of $a=0.13\leavevmode\nobreak\ $fm and with $\vec{q}=\vec{0}$.
We use 1000 eigenvectors of the even-odd Dirac operator (i.e. 2000
eigenvectors of the Dirac operator) to compute the low-mode contribution. The
results are depicted in Fig. 15.

We see that the split-even estimator outperforms the standard estimator, even
when the latter is combined with LMA. Less than 1000 inversions are needed to
reach the gauge noise. When the split-even estimator is combined with LMA,
only $O(10)$ sources are needed to reach the gauge noise. Combining LMA with
the split-even stochastic estimator, the computer time spent to compute those
disconnected loops is completely sub-dominant.

Figure 15: Square root of the variance of the vector loop function for the
four estimators described in the text. The lines represent a perfect scaling
of the variance assuming the variance is dominated by the stochastic noise.
The horizontal dashed line is the error from the low-mode contribution only
and is used as a normalization scale.

###  B.2 Contribution of the sub-leading diagrams to
$a_{\mu}^{\mathrm{HLbL};P}$

The second diagram in Fig. 2, with a single vector loop, is the only
disconnected contribution to the pion TFF. We emphasize that this
contribution, which vanishes in the SU(3) flavor limit, is included in the
main analysis. The goal of this section is to isolate its contribution to
$a_{\mu}^{\mathrm{HLbL};\pi^{0}}$. We start with the decomposition ${\cal
F}_{\pi^{0}\gamma^{*}\gamma^{*}}={\cal F}_{\pi^{0}\gamma^{*}\gamma^{*}}^{\rm
conn}+{\cal F}_{\pi^{0}\gamma^{*}\gamma^{*}}^{\rm disc}$ and treat the
disconnected contribution as a small correction. The result for ${\cal
F}_{\pi^{0}\gamma^{*}\gamma^{*}}^{\rm disc}$ on one ensemble is shown on left
panel of Fig. 16 where the data is fitted using the $z$-expansion described in
Section V.2 but setting $P=1$. We observe that this contribution decreases
much faster than the connected contribution at large virtualities. Finally,
the contribution to $a_{\mu}^{\mathrm{HLbL};\pi^{0}}$ is obtained from a
linearized version of the master formulae Eq. (16) and the results are shown
in the right panel of Fig. 16. In the continuum limit, and using a linear fit
in $a^{2}$ we find $-1.33(19)\times 10^{-11}$ compatible with the published
result $-1.0(0.5)\times 10^{-11}$ from [59], but with much smaller
uncertainties. Although this contribution is small, it cannot be neglected for
the pion at our level of precision.

Figure 16: Contribution of disconnected diagram with a single vector loop to
the pion transition form factor at our finest lattice spacing (left panel) and
to $a_{\mu}^{\mathrm{HLbL};\pi^{0}}$ in the continuum limit (right panel).

The last diagram in Fig. 2 contains a pseudoscalar loop and contributes only
to the $\eta$ and $\eta^{\prime}$ TFFs. The relative contribution is shown in
Fig. 4. For those states, the two contributions that involves a single vector
loop partly compensate each other and their contribution turns out to be
negligible at the current level of precision.

## References

  * Abi  _et al._ [2021] B. Abi _et al._ (Muon g-2), Measurement of the Positive Muon Anomalous Magnetic Moment to 0.46 ppm, Phys. Rev. Lett. 126, 141801 (2021), arXiv:2104.03281 [hep-ex] .
  * Bennett  _et al._ [2006] G. W. Bennett _et al._ (Muon g-2), Final Report of the Muon E821 Anomalous Magnetic Moment Measurement at BNL, Phys. Rev. D 73, 072003 (2006), arXiv:hep-ex/0602035 .
  * Aoyama  _et al._ [2020] T. Aoyama _et al._ , The anomalous magnetic moment of the muon in the Standard Model, Phys. Rept. 887, 1 (2020), arXiv:2006.04822 [hep-ph] .
  * Grange  _et al._ [2015] J. Grange _et al._ (Muon g-2), Muon (g-2) Technical Design Report, (2015), arXiv:1501.06858 [physics.ins-det] .
  * Abe  _et al._ [2019] M. Abe _et al._ , A New Approach for Measuring the Muon Anomalous Magnetic Moment and Electric Dipole Moment, PTEP 2019, 053C02 (2019), arXiv:1901.03047 [physics.ins-det] .
  * Della Morte  _et al._ [2017] M. Della Morte, A. Francis, V. Gülpers, G. Herdoíza, G. von Hippel, H. Horch, B. Jäger, H. B. Meyer, A. Nyffeler, and H. Wittig, The hadronic vacuum polarization contribution to the muon $g-2$ from lattice QCD, JHEP 10, 020, arXiv:1705.01775 [hep-lat] .
  * Chakraborty  _et al._ [2018] B. Chakraborty _et al._ (Fermilab Lattice, LATTICE-HPQCD, MILC), Strong-Isospin-Breaking Correction to the Muon Anomalous Magnetic Moment from Lattice QCD at the Physical Point, Phys. Rev. Lett. 120, 152001 (2018), arXiv:1710.11212 [hep-lat] .
  * Borsanyi  _et al._ [2018] S. Borsanyi _et al._ (Budapest-Marseille-Wuppertal), Hadronic vacuum polarization contribution to the anomalous magnetic moments of leptons from first principles, Phys. Rev. Lett. 121, 022002 (2018), arXiv:1711.04980 [hep-lat] .
  * 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), Calculation of the hadronic vacuum polarization contribution to the muon anomalous magnetic moment, Phys. Rev. Lett. 121, 022003 (2018), arXiv:1801.07224 [hep-lat] .
  * Giusti  _et al._ [2019a] D. Giusti, V. Lubicz, G. Martinelli, F. Sanfilippo, and S. Simula, Electromagnetic and strong isospin-breaking corrections to the muon $g-2$ from Lattice QCD+QED, Phys. Rev. D 99, 114502 (2019a), arXiv:1901.10462 [hep-lat] .
  * Shintani and Kuramashi [2019] E. Shintani and Y. Kuramashi (PACS), Hadronic vacuum polarization contribution to the muon $g-2$ with 2+1 flavor lattice QCD on a larger than (10 fm$)^{4}$ lattice at the physical point, Phys. Rev. D 100, 034517 (2019), arXiv:1902.00885 [hep-lat] .
  * Davies  _et al._ [2020] C. T. H. Davies _et al._ (Fermilab Lattice, LATTICE-HPQCD, MILC), Hadronic-vacuum-polarization contribution to the muon’s anomalous magnetic moment from four-flavor lattice QCD, Phys. Rev. D 101, 034512 (2020), arXiv:1902.04223 [hep-lat] .
  * Gérardin  _et al._ [2019a] A. Gérardin, M. Cè, G. von Hippel, B. Hörz, H. B. Meyer, D. Mohler, K. Ottnad, J. Wilhelm, and H. Wittig, The leading hadronic contribution to $(g-2)_{\mu}$ from lattice QCD with $N_{\rm f}=2+1$ flavours of O($a$) improved Wilson quarks, Phys. Rev. D 100, 014510 (2019a), arXiv:1904.03120 [hep-lat] .
  * Aubin  _et al._ [2020] C. Aubin, T. Blum, C. Tu, M. Golterman, C. Jung, and S. Peris, Light quark vacuum polarization at the physical point and contribution to the muon $g-2$, Phys. Rev. D 101, 014503 (2020), arXiv:1905.09307 [hep-lat] .
  * Giusti and Simula [2019] D. Giusti and S. Simula, Lepton anomalous magnetic moments in Lattice QCD+QED, PoS LATTICE2019, 104 (2019), arXiv:1910.03874 [hep-lat] .
  * Borsanyi  _et al._ [2021] S. Borsanyi _et al._ , Leading hadronic contribution to the muon magnetic moment from lattice QCD, Nature 593, 51 (2021), arXiv:2002.12347 [hep-lat] .
  * Lehner and Meyer [2020] C. Lehner and A. S. Meyer, Consistency of hadronic vacuum polarization between lattice QCD and the R-ratio, Phys. Rev. D 101, 074515 (2020), arXiv:2003.04177 [hep-lat] .
  * Giusti and Simula [2022] D. Giusti and S. Simula, Window contributions to the muon hadronic vacuum polarization with twisted-mass fermions, PoS LATTICE2021, 189 (2022), arXiv:2111.15329 [hep-lat] .
  * Wang  _et al._ [2022] G. Wang, T. Draper, K.-F. Liu, and Y.-B. Yang (chiQCD), Muon g-2 with overlap valence fermion, (2022), arXiv:2204.01280 [hep-lat] .
  * Aubin  _et al._ [2022] C. Aubin, T. Blum, M. Golterman, and S. Peris, Muon anomalous magnetic moment with staggered fermions: Is the lattice spacing small enough?, Phys. Rev. D 106, 054503 (2022), arXiv:2204.12256 [hep-lat] .
  * Davier  _et al._ [2017] M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, Reevaluation of the hadronic vacuum polarisation contributions to the Standard Model predictions of the muon $g-2$ and ${\alpha(m_{Z}^{2})}$ using newest hadronic cross-section data, Eur. Phys. J. C77, 827 (2017), arXiv:1706.09436 [hep-ph] .
  * Keshavarzi  _et al._ [2018] A. Keshavarzi, D. Nomura, and T. Teubner, Muon $g-2$ and $\alpha(M_{Z}^{2})$: a new data-based analysis, Phys. Rev. D97, 114025 (2018), arXiv:1802.02995 [hep-ph] .
  * Colangelo  _et al._ [2019] G. Colangelo, M. Hoferichter, and P. Stoffer, Two-pion contribution to hadronic vacuum polarization, JHEP 02, 006, arXiv:1810.00007 [hep-ph] .
  * Hoferichter  _et al._ [2019] M. Hoferichter, B.-L. Hoid, and B. Kubis, Three-pion contribution to hadronic vacuum polarization, JHEP 08, 137, arXiv:1907.01556 [hep-ph] .
  * Davier  _et al._ [2020] M. Davier, A. Hoecker, B. Malaescu, and Z. Zhang, A new evaluation of the hadronic vacuum polarisation contributions to the muon anomalous magnetic moment and to $\mathbf{\boldsymbol{\alpha}(m_{Z}^{2})}$, Eur. Phys. J. C80, 241 (2020), [Erratum: Eur. Phys. J. C80, 410 (2020)], arXiv:1908.00921 [hep-ph] .
  * Keshavarzi  _et al._ [2020a] A. Keshavarzi, D. Nomura, and T. Teubner, The $g-2$ of charged leptons, $\alpha(M_{Z}^{2})$ and the hyperfine splitting of muonium, Phys. Rev. D101, 014029 (2020a), arXiv:1911.00367 [hep-ph] .
  * Colangelo  _et al._ [2021] G. Colangelo, M. Hoferichter, and P. Stoffer, Constraints on the two-pion contribution to hadronic vacuum polarization, Phys. Lett. B 814, 136073 (2021), arXiv:2010.07943 [hep-ph] .
  * Crivellin  _et al._ [2020] A. Crivellin, M. Hoferichter, C. A. Manzari, and M. Montull, Hadronic Vacuum Polarization: $(g-2)_{\mu}$ versus Global Electroweak Fits, Phys. Rev. Lett. 125, 091801 (2020), arXiv:2003.04886 [hep-ph] .
  * Keshavarzi  _et al._ [2020b] A. Keshavarzi, W. J. Marciano, M. Passera, and A. Sirlin, Muon $g-2$ and $\Delta\alpha$ connection, Phys. Rev. D 102, 033002 (2020b), arXiv:2006.12666 [hep-ph] .
  * Malaescu and Schott [2021] B. Malaescu and M. Schott, Impact of correlations between $a_{\mu}$ and $\alpha_{\text{QED}}$ on the EW fit, Eur. Phys. J. C 81, 46 (2021), arXiv:2008.08107 [hep-ph] .
  * Blum  _et al._ [2015] T. Blum, S. Chowdhury, M. Hayakawa, and T. Izubuchi, Hadronic light-by-light scattering contribution to the muon anomalous magnetic moment from lattice QCD, Phys. Rev. Lett. 114, 012001 (2015), arXiv:1407.2923 [hep-lat] .
  * Blum  _et al._ [2016] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, and C. Lehner, Lattice Calculation of Hadronic Light-by-Light Contribution to the Muon Anomalous Magnetic Moment, Phys. Rev. D 93, 014503 (2016), arXiv:1510.07100 [hep-lat] .
  * Blum  _et al._ [2017] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, and C. Lehner, Connected and Leading Disconnected Hadronic Light-by-Light Contribution to the Muon Anomalous Magnetic Moment with a Physical Pion Mass, Phys. Rev. Lett. 118, 022005 (2017), arXiv:1610.04603 [hep-lat] .
  * Blum  _et al._ [2020] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, and C. Lehner, Hadronic Light-by-Light Scattering Contribution to the Muon Anomalous Magnetic Moment from Lattice QCD, Phys. Rev. Lett. 124, 132002 (2020), arXiv:1911.08123 [hep-lat] .
  * Blum  _et al._ [2023] T. Blum, N. Christ, M. Hayakawa, T. Izubuchi, L. Jin, C. Jung, C. Lehner, and C. Tu, Hadronic light-by-light contribution to the muon anomaly from lattice QCD with infinite volume QED at physical pion mass, (2023), arXiv:2304.04423 [hep-lat] .
  * Green  _et al._ [2015] J. Green, O. Gryniuk, G. von Hippel, H. B. Meyer, and V. Pascalutsa, Lattice QCD calculation of hadronic light-by-light scattering, Phys. Rev. Lett. 115, 222003 (2015), arXiv:1507.01577 [hep-lat] .
  * Green  _et al._ [2016] J. Green, N. Asmussen, O. Gryniuk, G. von Hippel, H. B. Meyer, A. Nyffeler, and V. Pascalutsa, Direct calculation of hadronic light-by-light scattering, PoS LATTICE2015, 109 (2016), arXiv:1510.08384 [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, Hadronic light-by-light scattering amplitudes from lattice QCD versus dispersive sum rules, Phys. Rev. D 98, 074501 (2018), arXiv:1712.00421 [hep-lat] .
  * Chao  _et al._ [2020] E.-H. Chao, A. Gérardin, J. R. Green, R. J. Hudspith, and H. B. Meyer, Hadronic light-by-light contribution to $(g-2)_{\mu}$ from lattice QCD with SU(3) flavor symmetry, 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, Hadronic light-by-light contribution to $(g-2)_{\mu}$ from lattice QCD: a complete calculation, Eur. Phys. J. C 81, 651 (2021), arXiv:2104.02632 [hep-lat] .
  * Asmussen  _et al._ [2022] N. Asmussen, E.-H. Chao, A. Gérardin, J. R. Green, R. J. Hudspith, H. B. Meyer, and A. Nyffeler, Hadronic light-by-light scattering contribution to the muon $g-2$ from lattice QCD: semi-analytical calculation of the QED kernel, (2022), arXiv:2210.12263 [hep-lat] .
  * Chao  _et al._ [2022a] E.-H. Chao, R. J. Hudspith, A. Gérardin, J. R. Green, and H. B. Meyer, The charm-quark contribution to light-by-light scattering in the muon $(g-2)$ from lattice QCD, Eur. Phys. J. C 82, 664 (2022a), arXiv:2204.08844 [hep-lat] .
  * Colangelo  _et al._ [2014a] G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer, Dispersive approach to hadronic light-by-light scattering, JHEP 09, 091, arXiv:1402.7081 [hep-ph] .
  * Colangelo  _et al._ [2014b] G. Colangelo, M. Hoferichter, B. Kubis, M. Procura, and P. Stoffer, Towards a data-driven analysis of hadronic light-by-light scattering, Phys. Lett. B 738, 6 (2014b), arXiv:1408.2517 [hep-ph] .
  * Colangelo  _et al._ [2015] G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer, Dispersion relation for hadronic light-by-light scattering: theoretical foundations, JHEP 09, 074, arXiv:1506.01386 [hep-ph] .
  * Colangelo  _et al._ [2017a] G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer, Rescattering effects in the hadronic-light-by-light contribution to the anomalous magnetic moment of the muon, Phys. Rev. Lett. 118, 232001 (2017a), arXiv:1701.06554 [hep-ph] .
  * Colangelo  _et al._ [2017b] G. Colangelo, M. Hoferichter, M. Procura, and P. Stoffer, Dispersion relation for hadronic light-by-light scattering: two-pion contributions, JHEP 04, 161, arXiv:1702.07347 [hep-ph] .
  * Chao  _et al._ [2022b] E.-H. Chao, R. J. Hudspith, A. Gérardin, J. R. Green, H. B. Meyer, and K. Ottnad, Hadronic light-by-light contribution to $(g-2)_{\mu}$ from lattice QCD: a complete calculation, PoS LATTICE2021, 209 (2022b).
  * Aubert  _et al._ [2009] B. Aubert _et al._ (BaBar), Measurement of the $\gamma\gamma^{*}\to\pi^{0}$ transition form factor, Phys. Rev. D 80, 052002 (2009), arXiv:0905.4778 [hep-ex] .
  * Behrend  _et al._ [1991] H. J. Behrend _et al._ (CELLO), A Measurement of the $\pi^{0}$, $\eta$ and $\eta^{\prime}$ electromagnetic form-factors, Z. Phys. C 49, 401 (1991).
  * Gronberg  _et al._ [1998] J. Gronberg _et al._ (CLEO), Measurements of the meson - photon transition form-factors of light pseudoscalar mesons at large momentum transfer, Phys. Rev. D 57, 33 (1998), arXiv:hep-ex/9707031 .
  * Uehara  _et al._ [2012] S. Uehara _et al._ (Belle), Measurement of $\gamma\gamma^{*}\to\pi^{0}$ transition form factor at Belle, Phys. Rev. D 86, 092007 (2012), arXiv:1205.3249 [hep-ex] .
  * del Amo Sanchez  _et al._ [2011] P. del Amo Sanchez _et al._ (BaBar), Measurement of the $\gamma\gamma^{*}\to\eta$ and $\gamma\gamma^{*}\to\eta^{\prime}$ transition form factors, Phys. Rev. D 84, 052001 (2011), arXiv:1101.1142 [hep-ex] .
  * Acciarri  _et al._ [1998] M. Acciarri _et al._ (L3), Measurement of eta-prime (958) formation in two photon collisions at LEP-1, Phys. Lett. B 418, 399 (1998).
  * Lees  _et al._ [2018] J. P. Lees _et al._ (BaBar), Measurement of the $\gamma^{\star}\gamma^{\star}\to\eta^{\prime}$ transition form factor, Phys. Rev. D 98, 112002 (2018), arXiv:1808.08038 [hep-ex] .
  * Tanabashi  _et al._ [2018] M. Tanabashi _et al._ (Particle Data Group), Review of Particle Physics, Phys. Rev. D 98, 030001 (2018).
  * Larin  _et al._ [2020] I. Larin _et al._ (PrimEx-II), Precision measurement of the neutral pion lifetime, Science 368, 506 (2020).
  * Gérardin  _et al._ [2016] A. Gérardin, H. B. Meyer, and A. Nyffeler, Lattice calculation of the pion transition form factor $\pi^{0}\to\gamma^{*}\gamma^{*}$, Phys. Rev. D 94, 074507 (2016), arXiv:1607.08174 [hep-lat] .
  * Gérardin  _et al._ [2019b] A. Gérardin, H. B. Meyer, and A. Nyffeler, Lattice calculation of the pion transition form factor with $N_{f}=2+1$ Wilson quarks, Phys. Rev. D 100, 034520 (2019b), arXiv:1903.09471 [hep-lat] .
  * Hoferichter  _et al._ [2018a] M. Hoferichter, B.-L. Hoid, B. Kubis, S. Leupold, and S. P. Schneider, Dispersion relation for hadronic light-by-light scattering: pion pole, JHEP 10, 141, arXiv:1808.04823 [hep-ph] .
  * Hoferichter  _et al._ [2018b] M. Hoferichter, B.-L. Hoid, B. Kubis, S. Leupold, and S. P. Schneider, Pion-pole contribution to hadronic light-by-light scattering in the anomalous magnetic moment of the muon, Phys. Rev. Lett. 121, 112002 (2018b), arXiv:1805.01471 [hep-ph] .
  * Alexandrou  _et al._ [2022] C. Alexandrou _et al._ , The $\eta\rightarrow\gamma^{*}\gamma^{*}$ transition form factor and the hadronic light-by-light $\eta$-pole contribution to the muon $g-2$ from lattice QCD, (2022), arXiv:2212.06704 [hep-lat] .
  * [63] BMW, Lattice QCD calculation of the $\eta$ and $\eta^{\prime}$ masses using staggered quarks, in preparation .
  * Ji and Jung [2001a] X.-d. Ji and C.-w. Jung, Studying hadronic structure of the photon in lattice QCD, Phys. Rev. Lett. 86, 208 (2001a), arXiv:hep-lat/0101014 .
  * Ji and Jung [2001b] X.-d. Ji and C.-w. Jung, Photon structure functions from quenched lattice QCD, Phys. Rev. D 64, 034506 (2001b), arXiv:hep-lat/0103007 .
  * Blossier  _et al._ [2009] B. Blossier, M. Della Morte, G. von Hippel, T. Mendes, and R. Sommer, On the generalized eigenvalue method for energies and matrix elements in lattice field theory, JHEP 04, 094, arXiv:0902.1265 [hep-lat] .
  * Jegerlehner and Nyffeler [2009] F. Jegerlehner and A. Nyffeler, The Muon g-2, Phys. Rept. 477, 1 (2009), arXiv:0902.3360 [hep-ph] .
  * Follana  _et al._ [2007] E. Follana, Q. Mason, C. Davies, K. Hornbostel, G. P. Lepage, J. Shigemitsu, H. Trottier, and K. Wong (HPQCD, UKQCD), Highly improved staggered quarks on the lattice, with applications to charm physics, Phys. Rev. D 75, 054502 (2007), arXiv:hep-lat/0610092 .
  * Altmeyer  _et al._ [1993] R. Altmeyer, K. D. Born, M. Gockeler, R. Horsley, E. Laermann, and G. Schierholz (MT(c)), The Hadron spectrum in QCD with dynamical staggered fermions, Nucl. Phys. B 389, 445 (1993).
  * Lepage [1998] P. Lepage, Perturbative improvement for lattice QCD: An Update, Nucl. Phys. B Proc. Suppl. 60, 267 (1998), arXiv:hep-lat/9707026 .
  * Lepage [1999] G. P. Lepage, Flavor symmetry restoration and Symanzik improvement for staggered quarks, Phys. Rev. D 59, 074502 (1999), arXiv:hep-lat/9809157 .
  * DeGrand and Schaefer [2005] T. A. DeGrand and S. Schaefer, Improving meson two-point functions by low-mode averaging, Nucl. Phys. B Proc. Suppl. 140, 296 (2005), arXiv:hep-lat/0409056 .
  * Giusti  _et al._ [2004] L. Giusti, P. Hernandez, M. Laine, P. Weisz, and H. Wittig, Low-energy couplings of QCD from current correlators near the chiral limit, JHEP 04, 013, arXiv:hep-lat/0402002 .
  * Bali  _et al._ [2010] G. S. Bali, S. Collins, and A. Schafer, Effective noise reduction techniques for disconnected loops in Lattice QCD, Comput. Phys. Commun. 181, 1570 (2010), arXiv:0910.3970 [hep-lat] .
  * Blum  _et al._ [2013] T. Blum, T. Izubuchi, and E. Shintani, New class of variance-reduction techniques using lattice symmetries, Phys. Rev. D 88, 094503 (2013), arXiv:1208.4349 [hep-lat] .
  * Shintani  _et al._ [2015] E. Shintani, R. Arthur, T. Blum, T. Izubuchi, C. Jung, and C. Lehner, Covariant approximation averaging, Phys. Rev. D 91, 114511 (2015), arXiv:1402.0244 [hep-lat] .
  * Gregory  _et al._ [2008] E. B. Gregory, A. C. Irving, C. M. Richards, and C. McNeile, Methods for Pseudoscalar Flavour-Singlet Mesons with Staggered Fermions, Phys. Rev. D 77, 065019 (2008), arXiv:0709.4224 [hep-lat] .
  * Giusti  _et al._ [2019b] L. Giusti, T. Harris, A. Nada, and S. Schaefer, Frequency-splitting estimators of single-propagator traces, Eur. Phys. J. C 79, 586 (2019b), arXiv:1903.10447 [hep-lat] .
  * Feng  _et al._ [2012] X. Feng, S. Aoki, H. Fukaya, S. Hashimoto, T. Kaneko, J.-i. Noaki, and E. Shintani, Two-photon decay of the neutral pion in lattice QCD, Phys. Rev. Lett. 109, 182001 (2012), arXiv:1206.1375 [hep-lat] .
  * Boyd  _et al._ [1996] C. G. Boyd, B. Grinstein, and R. F. Lebed, Model independent determinations of $\bar{B}\to D\ell\bar{\nu}$, $D^{*}\ell\bar{\nu}$ form-factors, Nucl. Phys. B 461, 493 (1996), arXiv:hep-ph/9508211 .
  * Lepage and Brodsky [1979] G. P. Lepage and S. J. Brodsky, Exclusive Processes in Quantum Chromodynamics: Evolution Equations for Hadronic Wave Functions and the Form-Factors of Mesons, Phys. Lett. B 87, 359 (1979).
  * Lepage and Brodsky [1980] G. P. Lepage and S. J. Brodsky, Exclusive Processes in Perturbative Quantum Chromodynamics, Phys. Rev. D 22, 2157 (1980).
  * Brodsky and Lepage [1981] S. J. Brodsky and G. P. Lepage, Large Angle Two Photon Exclusive Channels in Quantum Chromodynamics, Phys. Rev. D 24, 1808 (1981).
  * Nesterenko and Radyushkin [1983] V. A. Nesterenko and A. V. Radyushkin, Comparison of the QCD Sum Rule Approach and Perturbative QCD Analysis for $\gamma^{*}\gamma^{*}\to\pi^{0}$ Process, Sov. J. Nucl. Phys. 38, 284 (1983).
  * Novikov  _et al._ [1984] V. A. Novikov, M. A. Shifman, A. I. Vainshtein, M. B. Voloshin, and V. I. Zakharov, Use and Misuse of QCD Sum Rules, Factorization and Related Topics, Nucl. Phys. B 237, 525 (1984).
  * Ammer and Durr [2020] M. Ammer and S. Durr, Details of a staggered fermion data analysis, PoS LATTICE2019, 048 (2020), arXiv:1910.11046 [hep-lat] .
  * Bruno and Sommer [2023] M. Bruno and R. Sommer, On fits to correlated and auto-correlated data, Comput. Phys. Commun. 285, 108643 (2023), arXiv:2209.14188 [hep-lat] .
  * Masjuan and Sanchez-Puertas [2017] P. Masjuan and P. Sanchez-Puertas, Pseudoscalar-pole contribution to the $(g_{\mu}-2)$: a rational approach, Phys. Rev. D 95, 054026 (2017), arXiv:1701.05829 [hep-ph] .
  * Eichmann  _et al._ [2019] G. Eichmann, C. S. Fischer, E. Weil, and R. Williams, Single pseudoscalar meson pole and pion box contributions to the anomalous magnetic moment of the muon, Phys. Lett. B 797, 134855 (2019), [Erratum: Phys.Lett.B 799, 135029 (2019)], arXiv:1903.10844 [hep-ph] .
  * Larin  _et al._ [2011] I. Larin _et al._ (PrimEx), A New Measurement of the $\pi^{0}$ Radiative Decay Width, Phys. Rev. Lett. 106, 162303 (2011), arXiv:1009.1681 [nucl-ex] .
  * Workman  _et al._ [2022] R. L. Workman _et al._ (Particle Data Group), Review of Particle Physics, PTEP 2022, 083C01 (2022).
  * Browman  _et al._ [1974] A. Browman, J. DeWire, B. Gittelman, K. M. Hanson, E. Loh, and R. Lewis, The Radiative Width of the eta Meson, Phys. Rev. Lett. 32, 1067 (1974).
  * Bartel  _et al._ [1985] W. Bartel _et al._ (JADE), A measurement of the $\eta$ radiative width $\Gamma_{\eta\to\gamma\gamma}$, Phys. Lett. B 158, 511 (1985).
  * Williams  _et al._ [1988] D. Williams _et al._ (Crystal Ball), Formation of the Pseudoscalars $\pi^{0}$, $\eta$ and $\eta^{\prime}$ in the Reaction $\gamma\gamma\to\gamma\gamma$, Phys. Rev. D 38, 1365 (1988).
  * Roe  _et al._ [1990] N. A. Roe _et al._ , A Measurement of the Radiative Width of the $\eta$ and $\eta^{\prime}$ Mesons With the Asp Detector, Phys. Rev. D 41, 17 (1990).
  * Baru  _et al._ [1990] S. E. Baru _et al._ , Measurement of two photon widths of the A2, $\eta^{\prime}$, $\eta$, Z. Phys. C 48, 581 (1990).
  * Babusci  _et al._ [2013] D. Babusci _et al._ (KLOE-2), Measurement of $\eta$ meson production in $\gamma\gamma$ interactions and $\Gamma(\eta\to\gamma\gamma)$ with the KLOE detector, JHEP 01, 119, arXiv:1211.1845 [hep-ex] .
  * Liping Gan [2022] Liping Gan, The $\eta$ and $\eta^{\prime}$ physics at JLab, in  _Fifth Plenary Workshop of the Muon $g-2$ Theory Initiative_ (2022).
  * Nyffeler [2016] A. Nyffeler, Precision of a data-driven estimate of hadronic light-by-light scattering in the muon $g-2$: Pseudoscalar-pole contribution, Phys. Rev. D 94, 053006 (2016), arXiv:1602.03398 [hep-ph] .
  * Raya  _et al._ [2020] K. Raya, A. Bashir, and P. Roig, Contribution of neutral pseudoscalar mesons to $a_{\mu}^{HLbL}$ within a Schwinger-Dyson equations approach to QCD, Phys. Rev. D 101, 074021 (2020), arXiv:1910.05960 [hep-ph] .
  * Meng  _et al._ [2021] Y. Meng, X. Feng, C. Liu, T. Wang, and Z. Zou, First-principle calculation of $\eta_{c}\rightarrow 2\gamma$ decay width from lattice QCD, (2021), arXiv:2109.09381 [hep-lat] .
  * Feng and Jin [2019] X. Feng and L. Jin, QED self energies from lattice QCD without power-law finite-volume errors, Phys. Rev. D 100, 094509 (2019), arXiv:1812.09817 [hep-lat] .
  * Gülpers  _et al._ [2014] V. Gülpers, A. Francis, B. Jäger, H. Meyer, G. von Hippel, and H. Wittig, The leading disconnected contribution to the anomalous magnetic moment of the muon, PoS LATTICE2014, 128 (2014), arXiv:1411.7592 [hep-lat] .