Lattice calculation of the pion transition form factor with $N_f=2+1$ Wilson quarks

TL;DR

This paper presents a lattice QCD calculation of the pion transition form factor with $N_f=2+1$ Wilson quarks, improving accuracy for muon g-2 hadronic light-by-light scattering estimates by including strange quark effects and reducing discretization errors.

Contribution

It introduces a more precise, model-independent lattice estimate of the pion-pole contribution to muon g-2, incorporating improvements like on-shell $ ext{O}(a)$-improvement and extended kinematic coverage.

Findings

Normalized the form factor with 3.5% precision.

Confirmed low-energy constant estimates in chiral perturbation theory.

Provided a lattice estimate of the pion-pole contribution to muon g-2: (59.7 ± 3.6) × 10^{-11}.

Abstract

We present a lattice QCD calculation of the double-virtual neutral pion transition form factor, with the goal to cover the kinematic range relevant to hadronic light-by-light scattering in the muon $g-2$. Several improvements have been made compared to our previous work. First, we take into account the effects of the strange quark by using the $N_f=2+1$ CLS gauge ensembles. Secondly, we have implemented the on-shell $\mathcal{O}(a)$-improvement of the vector current to reduce the discretization effects associated with Wilson quarks. Finally, in order to have access to a wider range of photon virtualities, we have computed the transition form factor in a moving frame as well as in the pion rest-frame. After extrapolating the form factor to the continuum and to physical quark masses, we compare our results with phenomenology. We extract the normalization of the form factor with a precision of 3.5\% and confirm within our uncertainty previous somewhat conflicting estimates for a low-energy constant that appears in chiral perturbation theory for the decay $\pi^0 \to \gamma\gamma$ at NLO. With additional input from experiment and theory, we reproduce recent estimates for the decay width $\Gamma(\pi^0 \to \gamma\gamma)$. We also study the asymptotic large-$Q^2$ behavior of the transition form factor in the double-virtual case. Finally, we provide as our main result a more precise model-independent lattice estimate of the pion-pole contribution to hadronic light-by-light scattering in the muon $g-2$: $a_{\mu}^{\mathrm{HLbL}; \pi^0} = (59.7 \pm 3.6) \times 10^{-11}$. Using in addition the normalization of the form factor obtained by the PrimEx experiment, we get the lattice and data-driven estimate $a_{\mu}^{\mathrm{HLbL}; \pi^0} = (62.3 \pm 2.3) \times 10^{-11}$.

Tables

Table 1. Table 1: Parameters of the simulations: the bare coupling β = 6 / g 0 2 𝛽 6 superscript subscript 𝑔 0 2 \beta=6/g_{0}^{2} , the lattice resolution, the lattice spacing a 𝑎 a in physical units extracted from Bruno:2016plf , the light and strange hopping parameters κ l subscript 𝜅 𝑙 \kappa_{l} and κ s subscript 𝜅 𝑠 \kappa_{s} , the pion mass m π subscript 𝑚 𝜋 m_{\pi} , the ground state vector mass, the number of gauge configurations and the number of sources per configuration for the three-point correlation function. Ensembles with an asterisk are not included in the final analysis but used to control finite size-effects.

Id	$\beta$	$L^3\times T$	$a\left[\mathrm{fm}\right]$	${\kappa }_l$	${\kappa }_s$	$m_{\pi }\left[\mathrm{MeV}\right]$	$m_V\left[\mathrm{MeV}\right]$	$m_{\pi }L$	$\#$confs	$\#$src
H101	$3.40$	${32}^3\times 96$	0.08636	0.136760	0.13675962	416(6)	844(11)	5.8	1000	10
H102		${32}^3\times 96$		0.136865	0.13654934	354(5)	832(09)	5.0	1900	10
H105∗		${32}^3\times 96$		0.136970	0.13634079	281(4)	759(18)	3.9	2800	10
N101		${48}^3\times 128$		0.136970	0.13634079	280(4)	774(08)	5.9	1600	10
C101		${48}^3\times 96$		0.137030	0.13622204	224(3)	741(10)	4.7	2200	15
S400	$3.46$	${32}^3\times 128$	0.07634	0.136984	0.13670239	349(5)	821(09)	4.3	1700	20
N401		${48}^3\times 128$		0.137062	0.13654808	286(4)	793(09)	5.3	950	10
H200∗	$3.55$	${32}^3\times 96$	0.06426	0.137000	0.137000	419(6)	875(17)	4.4	2000	10
N202		${48}^3\times 128$		0.137000	0.137000	411(5)	859(10)	6.4	900	5
N203		${48}^3\times 128$		0.137080	0.13684028	346(5)	830(09)	5.4	1500	10
N200		${48}^3\times 128$		0.137140	0.13672086	284(3)	805(13)	4.4	1700	10
D200		${64}^3\times 128$		0.137200	0.13660175	200(3)	740(14)	4.2	1100	30
N300	$3.70$	${48}^3\times 128$	0.04981	0.137000	0.137000	422(5)	897(12)	5.1	1200	5
N302		${48}^3\times 128$		0.137064	0.13687218	343(5)	856(16)	4.2	1100	10
J303		${64}^3\times 192$		0.137123	0.13675466	258(3)	796(09)	4.2	650	10

Table 2. Table 2: Ground state energies E π ​ ( p → ) subscript 𝐸 𝜋 → 𝑝 E_{\pi}(\vec{p}) for the two values of the pion momentum and overlap factors Z π subscript 𝑍 𝜋 Z_{\pi} for each lattice ensemble. We also give the PCAC mass and the pion decay constant in lattice units without finite-size effect correction.

Id	$i{Z}_{\pi }$	$a{E}_{\pi }\left(\overrightarrow{0}\right)$	$a{E}_{\pi }\left(\overrightarrow{p}\right)$	$a{m}_{\mathrm{PCAC}}$	$a\sqrt{2}{f}_{\pi }$
H101	$0.2139\left(13\right)$	$0.1821\left(07\right)$	$0.2671\left(07\right)$	$0.009177\left(45\right)$	$0.06458\left(29\right)$
H102	$0.2086\left(19\right)$	$0.1550\left(09\right)$	$0.2492\left(09\right)$	$0.006499\left(51\right)$	$0.06151\left(27\right)$
H105	$0.2004\left(23\right)$	$0.1230\left(12\right)$	$0.2316\left(11\right)$	$0.003985\left(57\right)$	$0.05802\left(39\right)$
N101	$0.2011\left(15\right)$	$0.1224\left(05\right)$	$0.1792\left(05\right)$	$0.003990\left(32\right)$	$0.05832\left(30\right)$
C101	$0.1947\left(11\right)$	$0.0982\left(06\right)$	$0.1629\left(05\right)$	$0.002435\left(29\right)$	$0.05534\left(37\right)$
S400	$0.1606\left(16\right)$	$0.1352\left(08\right)$	$0.2387\left(09\right)$	$0.005684\left(27\right)$	$0.05463\left(21\right)$
N401	$0.1578\left(10\right)$	$0.1106\left(05\right)$	$0.1714\left(05\right)$	$0.003770\left(28\right)$	$0.05327\left(17\right)$
H200	$0.1160\left(12\right)$	$0.1365\left(07\right)$	$0.2382\left(10\right)$	$0.006863\left(24\right)$	$0.04805\left(27\right)$
N202	$0.1155\left(10\right)$	$0.1337\left(05\right)$	$0.1870\left(06\right)$	$0.006866\left(15\right)$	$0.04884\left(18\right)$
N203	$0.1135\left(08\right)$	$0.1127\left(04\right)$	$0.1728\left(05\right)$	$0.004749\left(17\right)$	$0.04699\left(16\right)$
N200	$0.1094\left(08\right)$	$0.0924\left(05\right)$	$0.1599\left(05\right)$	$0.003145\left(15\right)$	$0.04454\left(18\right)$
D200	$0.1044\left(08\right)$	$0.0650\left(04\right)$	$0.1173\left(05\right)$	$0.001554\left(12\right)$	$0.04259\left(17\right)$
N300	$0.0700\left(07\right)$	$0.1066\left(04\right)$	$0.1685\left(06\right)$	$0.005507\left(09\right)$	$0.03811\left(14\right)$
N302	$0.0646\left(08\right)$	$0.0867\left(05\right)$	$0.1556\left(09\right)$	$0.003725\left(10\right)$	$0.03570\left(22\right)$
J303	$0.0632\left(05\right)$	$0.0652\left(02\right)$	$0.1177\left(04\right)$	$0.002056\left(07\right)$	$0.03412\left(15\right)$

Table 3. Table 3: Coefficients of the z 𝑧 z -expansion, in GeV − 1 superscript GeV 1 \mathrm{GeV}^{-1} , defined through Eq. ( 43 ) in the continuum, at physical quark masses, for different values of N 𝑁 N . We use Q max 2 = 4 ​ GeV 2 superscript subscript 𝑄 2 4 superscript GeV 2 Q_{\max}^{2}=4\leavevmode\nobreak\ \mathrm{GeV}^{2} . The chi-squares per degree of freedom are respectively χ 2 / d . o . f . = 1.5 ,  1.2 ,  1.1 \chi^{2}/\mathrm{d.o.f.}=1.5,\,1.2,\,1.1 and for N = 3 𝑁 3 N=3 , the associated correlation matrix is given by Eq. ( 45 ).

$N$	$c_{00}$	$c_{01}$	$c_{11}$	$c_{20}$	$c_{21}$	$c_{22}$	$c_{30}$	$c_{31}$	$c_{32}$	$c_{33}$
1	0.2346(65)	$-0.0590\left(39\right)$	0.074(19)
2	0.2350(61)	$-0.0651\left(49\right)$	$-0.284\left(68\right)$	0.106(33)	0.109(46)	$-0.29\left(12\right)$
3	0.2345(63)	$-0.0746\left(52\right)$	$-0.338\left(86\right)$	0.145(43)	0.008(127)	$-0.92\left(55\right)$	0.34(10)	0.25(27)	$-1.27\left(79\right)$	1.16(1.40)

Table 4. Table 4: Normalization of the transition form factor obtained from the two sequences of Canterbury approximants C N N subscript superscript 𝐶 𝑁 𝑁 C^{N}_{N} and C N + 1 N subscript superscript 𝐶 𝑁 𝑁 1 C^{N}_{N+1} . As before, chi-square corresponds to uncorrelated fits.

	$C_1^0$	$C_2^1$	$C_1^1$
$\alpha \left[{\mathrm{GeV}}^{-1}\right]$	0.325(8)	0.253(12)	0.268(10)
${\chi }^2/d.o.f.$	4.1	1.2	1.3

Table 5. Table 5: Value of the higher twist coefficient δ 2 superscript 𝛿 2 \delta^{2} as a function of Q min 2 subscript superscript 𝑄 2 min Q^{2}_{\rm min} .

$Q_{\min }^2\left[{\mathrm{GeV}}^2\right]$	1.5	1.75	2.0	2.25
${\delta }^2\left[{\mathrm{GeV}}^2\right]$	0.13(5)	0.14(6)	0.15(6)	0.15(7)

Table 6. Table 6: Results of the fits of the LMD+V model using the modified double z 𝑧 z -expansion. The model is defined using Eq. ( 48c ) and the fit parameters in Eq. ( 53 ). The exact values are α = 0.264 ​ GeV − 1 𝛼 0.264 superscript GeV 1 \alpha=0.264\leavevmode\nobreak\ \mathrm{GeV}^{-1} , b π = 1.62 ​ GeV − 2 subscript 𝑏 𝜋 1.62 superscript GeV 2 b_{\pi}=1.62\leavevmode\nobreak\ \mathrm{GeV}^{-2} and a μ HLbL ; π 0 = 59.2 × 10 − 11 superscript subscript 𝑎 𝜇 HLbL superscript 𝜋 0 59.2 superscript 10 11 a_{\mu}^{\mathrm{HLbL};\pi^{0}}=59.2\times 10^{-11} . Results for a μ HLbL ; π 0 superscript subscript 𝑎 𝜇 HLbL superscript 𝜋 0 a_{\mu}^{\mathrm{HLbL};\pi^{0}} are given in units of 10 − 11 superscript 10 11 10^{-11} . The last column, d max subscript 𝑑 max d_{\rm max} , corresponds to the maximum deviation in percent between the exact TFF and the fit.

	$Q_{\max }^2=1{\mathrm{GeV}}^2$			$Q_{\max }^2=2{\mathrm{GeV}}^2$			$Q_{\max }^2=4{\mathrm{GeV}}^2$
$N$	$\alpha \left[{\mathrm{GeV}}^{-1}\right]$	$b_{\pi }\left[{\mathrm{GeV}}^{-2}\right]$	$a_{\mu }^{\mathrm{HLbL};{\pi }^0}$	$\alpha \left[{\mathrm{GeV}}^{-1}\right]$	$b_{\pi }\left[{\mathrm{GeV}}^{-2}\right]$	$a_{\mu }^{\mathrm{HLbL};{\pi }^0}$	$\alpha \left[{\mathrm{GeV}}^{-1}\right]$	$b_{\pi }\left[{\mathrm{GeV}}^{-2}\right]$	$a_{\mu }^{\mathrm{HLbL};{\pi }^0}$	$d_{\max }$
$1$	0.260	1.51	55.4	0.253	1.51	54.4	0.232	1.43	49.1	12.2%
$2$	0.265	1.63	57.0	0.264	1.55	58.2	0.262	1.40	58.4	3.5%
$3$	0.264	1.62	57.9	0.264	1.65	57.9	0.265	1.71	58.3	0.6%
$4$	0.264	1.62	57.8	0.264	1.61	58.3	0.264	1.60	59.2	0.2%
$5$	0.264	1.62	58.4	0.264	1.63	59.4	0.264	1.62	59.8	0.06%

Table 7. Table 7: Coefficients of the z 𝑧 z -expansion, in GeV − 1 superscript GeV 1 \mathrm{GeV}^{-1} , defined in Eq. ( 43 ) and obtained from a fit on a single ensemble using local vector currents at the source and at the sink. The vector meson masses used in Eq. ( 41 ) are given in Table  1 .

Id	$c_{00}$	$c_{01}$	$c_{11}$	$c_{20}$	$c_{21}$
N101	0.2456(46)	$-0.0755\left(61\right)$	$-0.265\left(99\right)$	0.065(54)	0.02(14)
N200	0.2301(44)	$-0.0675\left(43\right)$	$-0.280\left(64\right)$	0.071(32)	0.09(11)
D200	0.2484(48)	$-0.0728\left(51\right)$	$-0.399\left(94\right)$	0.148(47)	$-0.02\left(14\right)$

Equations

$$M_{\mu \nu }(p,q_1)=i\int {\mathrm{d}}^{4}x{e}^{i{q}_1\cdot x}⟨\Omega |T\{J_{\mu }(x)J_{\nu }(0)\}|{\pi }^0(p)⟩={ϵ}_{\mu \nu \alpha \beta }{q}_1^{\alpha }{q}_2^{\beta }{\mathcal{F}}_{{\pi }^0{\gamma }^{\ast }{\gamma }^{\ast }}(q_1^2,q_2^2),$$

$$M_{\mu \nu }=(i^{n_0})M_{\mu \nu }^{\mathrm{E}},M_{\mu \nu }^{\mathrm{E}}\equiv -\int \mathrm{d}\tau e^{{\omega }_1\tau }\int {\mathrm{d}}^{3}z{e}^{-i{\vec{q}}_1\vec{z}}⟨0|T\left\{J_{\mu }(\vec{z},\tau )J_{\nu }(\vec{0},0)\right\}|\pi (p)⟩,$$

$$M_{\mu \nu }^{\mathrm{E}}=\frac{2E_{\pi }}{Z_{\pi }}{\int }_{-\infty }^{\infty }\mathrm{d}\tau e^{{\omega }_1\tau }{\tilde{A}}_{\mu \nu }(\tau ),$$

$${\tilde{A}}_{\mu \nu }(\tau )\equiv \underset{t_{\pi }\to +\infty }{\lim }{e}^{E_{\pi }(t_f-t_0)}{C}_{\mu \nu }^{(3)}(\tau ,t_{\pi }).$$

$$C^{(3)}_{\mu\nu}(\tau,t_{\pi})\equiv a^{6}\sum_{\vec{x},\vec{z}}\,\big{\langle}J_{\mu}(\vec{z},t_{i})J_{\nu}(\vec{0},t_{f})P^{{\dagger}}(\vec{x},t_{0})\big{\rangle}\,e^{i\vec{p}\,\vec{x}}\,e^{-i\vec{q}_{1}\vec{z}}\,,$$

$$\begin{array}{rl}{q}_1^2& ={\omega }_1^2-{\vec{q}}_1^2\\ q_2^2& =(E_{\pi }-{\omega }_1{)}^2-(\vec{p}-{\vec{q}}_1{)}^2\end{array}\mathrm{with}{\vec{q}}_1=\frac{2\pi }{L}\vec{n},\vec{n}\in {\mathbb{Z}}^3.$$

$$M_{\mu \nu }^{\mathrm{E}}=\frac{2E_{\pi }}{Z_{\pi }}{\int }_{-\infty }^{\infty }\mathrm{d}\tau e^{i\tilde{\omega }\tau }{\tilde{A}}_{\mu \nu }(\tau ).$$

$${\tilde{A}}_{\mu \nu }(\tau )=\frac{Z_{\pi }}{4\pi E_{\pi }}{\int }_{-\infty }^{\infty }\mathrm{d}\tilde{\omega }{M}_{\mu \nu }^{\mathrm{E}}{e}^{-i\tilde{\omega }\tau }.$$

$$q_{\mu \nu }\equiv {ϵ}_{\mu \nu \alpha \beta }{q}_1^{\alpha }{q}_2^{\beta }=P_{\mu \nu }{\omega }_1+Q_{\mu \nu },$$

$${\tilde{A}}_{\mu \nu }(\tau )=-i{Q}_{\mu \nu }^E{\tilde{A}}^{(1)}(\tau )+P_{\mu \nu }^E\frac{\mathrm{d}{\tilde{A}}^{(1)}}{\mathrm{d}\tau }(\tau ),$$

$${\tilde{A}}^{(1)}(\tau )=\frac{i{Z}_{\pi }}{4\pi E_{\pi }}{\int }_{-\infty }^{\infty }\mathrm{d}\tilde{\omega }{\mathcal{F}}_{{\pi }^0{\gamma }^{\ast }{\gamma }^{\ast }}(q_1^2,q_2^2)e^{-i\tilde{w}\tau }.$$

$${\tilde{A}}_{0k}(\tau )$$

$${ϵ}^{\prime k}{\tilde{A}}_{kl}(\tau ){ϵ}^l$$

$${\tilde{A}}_{12}(\tau )\equiv -i{E}_{\pi }{p}_z{\tilde{A}}^{(2)}(\tau ),\vec{p}=(2\pi /L)\hat{z}.$$

$${\tilde{A}}_{\mu \nu }(\tau ;{\vec{q}}_1,{\vec{q}}_2)={\tilde{A}}_{\nu \mu }(-\tau ;{\vec{q}}_2,{\vec{q}}_1)e^{-E_{\pi }\tau },$$

$$V_{\mu }^{l,a}(x)$$

$$V_{\mu }^{c,a}(x)$$

$$V_{\mu }^I(x)=V_{\mu }(x)+a{c}_{\mathrm{V}}(g_0){\partial }_{\nu }{\Sigma }_{\mu \nu }(x)$$

$$\mathrm{tr}(\lambda V_{\mu})_{R}=Z_{\mathrm{V}}(\widetilde{g}_{0})\Big{[}\left(1+3\,\overline{b}_{\rm V}(g_{0})\,am_{\rm q}^{\rm av}\right)\mathrm{tr}(\lambda V_{\mu}^{I})+\frac{1}{2}b_{\rm V}(g_{0})\,\mathrm{tr}(\{\lambda,aM_{\rm q}\}V_{\mu}^{I})\\ +f_{\rm V}(g_{0})\,\mathrm{tr}(\lambda\,aM_{\rm q})\,\mathrm{tr}(V_{\mu}^{I})\Big{]}$$

$$J_{\mu }={\hat{V}}_{\mu }^3+\frac{1}{\sqrt{3}}{\hat{V}}_{\mu }^8,$$

$${\hat{V}}_{\mu }^3$$

$${\hat{V}}_{\mu }^8$$

$$Z_3$$

$$Z_8$$

$$Z_{80}$$

$$⟨J_{\mu }^l(z)J_{\nu }^l(y)P^{†}(x)⟩=\frac{\sqrt{2}}{3}{Z}_3\left(Z_8+\sqrt{3}{Z}_{80}\right)\mathrm{Re}\mathrm{Tr}\left[G_l(x,z){\gamma }_{\mu }{G}_l(z,y){\gamma }_{\nu }{G}_l(y,x){\gamma }_5\right]$$

$$-\frac{\sqrt{2}}{6}{Z}_3\left[\left(Z_8+\sqrt{3}{Z}_{80}\right)\mathrm{Tr}\left[G_l(y,y){\gamma }_{\nu }\right]-\left(Z_8-\frac{\sqrt{3}}{2}{Z}_{80}\right)\mathrm{Tr}\left[G_s(y,y){\gamma }_{\nu }\right]\right]\mathrm{Tr}\left[G_l(z,x){\gamma }_5{G}_l(x,z){\gamma }_{\mu }\right]$$

$$-\frac{\sqrt{2}}{6}{Z}_3\left[\left(Z_8+\sqrt{3}{Z}_{80}\right)\mathrm{Tr}\left[G_l(z,z){\gamma }_{\mu }\right]-\left(Z_8-\frac{\sqrt{3}}{2}{Z}_{80}\right)\mathrm{Tr}\left[G_s(z,z){\gamma }_{\mu }\right]\right]\mathrm{Tr}\left[G_l(y,x){\gamma }_5{G}_l(x,y){\gamma }_{\nu }\right],$$

$$⟨J_{\mu }^l(z)J_{\nu }^c(y)P^{†}(x)⟩=\frac{\sqrt{2}}{6}\left(Z_3+Z_8+\sqrt{3}{Z}_{80}\right)\mathrm{Re}\mathrm{Tr}\left[G_l(x,z){\gamma }_{\mu }{G}_l(z,y^{\prime }){\mathcal{V}}_{\nu }^{(y)}(y^{\prime },y^{\prime \prime })G_l(y^{\prime \prime },x){\gamma }_5\right]$$

$$-\frac{\sqrt{2}}{6}\left[\left(Z_8+\sqrt{3}{Z}_{80}\right)\mathrm{Tr}\left[G_l(z,z){\gamma }_{\mu }\right]-\left(Z_8-\frac{\sqrt{3}}{2}{Z}_{80}\right)\mathrm{Tr}\left[G_s(z,z){\gamma }_{\mu }\right]\right]$$

$$\times \mathrm{Tr}\left[G_l(y^{\prime \prime },x){\gamma }_5{G}_l(x,y^{\prime }){\mathcal{V}}_{\nu }^{(y)}(y^{\prime },y^{\prime \prime })\right]$$

$$\displaystyle-\frac{\sqrt{2}}{6}Z_{3}\Big{[}\mathrm{Tr}\left[G_{l}(y^{\prime\prime},y^{\prime})\mathcal{V}^{(y)}_{\nu}(y^{\prime},y^{\prime\prime})\right]-\mathrm{Tr}\left[G_{s}(y^{\prime\prime},y^{\prime})\mathcal{V}^{(y)}_{\nu}(y^{\prime},y^{\prime\prime})\right]\Big{]}\mathrm{Tr}\left[G_{l}(z,x)\gamma_{5}G_{l}(x,z)\gamma_{\mu}\right]\,,$$

$${\mathcal{V}}_{\mu }^{(y)}(y^{\prime },y^{\prime \prime })G(y^{\prime \prime },x)=\frac{1}{2}\left({\delta }_{y^{\prime },y+a\hat{\mu }}[1+{\gamma }_{\mu }]U_{\mu }^{†}(y)G(y,x)-{\delta }_{y^{\prime },y}[1-{\gamma }_{\mu }]U_{\mu }(y)G(y+a\hat{\mu },x)\right).$$

$${\partial }_{\mu }^{s}f(x)=\frac{f(x+a)-f(x-a)}{2a},{\partial }_{\mu }^{b}f(x)=\frac{f(x)-f(x-a)}{a},{\partial }_{\mu }^{f}f(x)=\frac{f(x+a)-f(x)}{a}.$$

$$\Delta {\tilde{A}}^{(1)}(\tau )=\frac{|{\tilde{A}}^{(1),\mathrm{lc}}(\tau )-{\tilde{A}}^{(1),\mathrm{ll}}(\tau )|}{\delta {\tilde{A}}^{(1),\mathrm{lc}}(\tau )},$$

Full text

††preprint: MITP/19-014

I Introduction

There is a long-standing discrepancy between the Standard Model estimate of the muon anomalous magnetic moment and its experimental determination pdg:2018 . Two new experiments, E989 at Fermilab Venanzoni:2014ixa and E34 at J-PARC Otani:2015lra , plan to reduce the experimental error by a factor of 4 in the near future. The theory error is completely dominated by hadronic contributions: the hadronic vacuum polarization (HVP), which enters at order $\alpha_{\rm e}^{2}$ in the fine-structure constant $\alpha_{\rm e}$, and the hadronic light-by-light (HLbL) scattering at order $\alpha_{\rm e}^{3}$. The former is usually obtained using dispersive methods which rely on the e^{+}e^{-}\to\leavevmode\nobreak\hadrons cross sections, accessible from experiments Davier:2017zfy ; Anastasi:2017eio . Thus, a more accurate determination essentially relies on the availability of precise measurements. Lattice QCD is also a promising tool, and has made a lot of progress in recent years. Even if not yet competitive with the dispersive approach, it became a mature field where most sources of systematic uncertainties have now been addressed, see e.g. Ref. Meyer:2018til for a recent review. As for the HLbL contribution, the situation is less favorable. Until recently, all estimates for the HLbL contribution were based on model calculations, leading to the Glasgow consensus Prades:2009tw (see also Jegerlehner:2009ry ), for which errors are difficult to estimate. The single largest contribution to $a_{\mu}^{\rm HLbL}$ is given by the pion-pole contribution Knecht:2001qf with a prescription for its evaluation that has been confirmed in the recently proposed dispersive approaches to HLbL HLbL_DR ; Pauk:2014jza . Its determination relies on the knowledge of the neutral pion transition form factor (TFF). Two groups have also started the direct calculation of the HLbL scattering contribution, from first principles, using lattice QCD Blum:2014oka ; Blum:2015gfa ; Blum:2016lnc ; Green:2015mva ; Asmussen:2016lse ; Asmussen:2017bup . In the Mainz approach, the calculation involves the convolution of a QED kernel function, computed semianalytically in infinite, continuous coordinate space, and a QCD four-point correlation function, computed on the lattice. The long-distance contribution is expected to be dominated by the pion-pole contribution. However, this region suffers from large statistical errors and is also more affected by finite-size effects. The TFF is therefore a key ingredient to first reduce the statistical error by constraining the tail of the integrand at large distances and, second, to estimate and correct for the dominant finite-size effects due to pions.

The TFF is also interesting from a theoretical point of view. First, in the low energy region, it is directly related to the Adler-Bell-Jackiw (ABJ) chiral anomaly Adler:1969gk ; Bell:1969ts . Secondly, for large virtualities, the asymptotic behavior of the TFF is predicted by the Brodsky-Lepage analysis and the operator product expansion (OPE) in the single- and double-virtual case respectively BL_3_papers ; Nesterenko:1982dn ; Novikov:1983jt . Testing these predictions is a remarkable test of QCD over a large range of length scales.

The normalization of the TFF has been measured by the PrimEx experiment with a precision of 1.4 % Larin:2010kq and this error should be reduced by a factor of 2 by the PrimeEx-II experiment Gasparian:2016oyl ; Primex2Mainz . At finite virtualities, experimental data are only available in the single-virtual case and for rather large virtualities above $0.6\leavevmode\nobreak\ \mathrm{GeV}^{2}$ exp but the BESIII experiment plans to have data in the region 0.3 - 3 $\mathrm{GeV}^{2}$, relevant for the muon $g-2$, in the near future Denig:2014mma . Finally, no experimental data exist to date in the regime of two virtual photons.

In our previous work Gerardin:2016cqj , we have shown that lattice QCD can provide a precise estimate of the TFF in the full spacelike region relevant for the pion-pole contribution to the muon $(g-2)$. This work is an update of our previous lattice calculation of the neutral pion TFF Gerardin:2016cqj , and includes several major improvements. First, it is based on gauge configurations with $N_{f}=2+1$ dynamical flavors, which have been generated as part of the CLS initiative. Second, on-shell $\mathcal{O}(a)$-improvement of correlation functions has been implemented to reduce discretization effects when approaching the continuum limit. Finally, in addition to the pion rest frame, a new frame where the pion carries one unit of momentum, typically in the range of 300 to 400 MeV, is considered. This allows us to probe larger photon virtualities, especially in the single-virtual case, where virtualities as high as $1.5\leavevmode\nobreak\ \mathrm{GeV}^{2}$ are now available. In all cases, we probe a much wider range of virtualities than in our previous study and have increased statistics significantly. We also address many sources of systematic errors, including finite-size effects, hypercubic lattice artifacts due to the broken rotational invariance down to the isometry group H(3) on the lattice and disconnected contributions.

As a benchmark of our calculation, we reproduce the anomaly constraint in the continuum and at the physical pion mass with a precision of 3.5%. From the pion mass dependence of the normalization of the form factor we extract the corresponding low-energy constant (LEC) appearing in the chiral Lagrangian needed for the next-to-leading order (NLO) calculation of the decay $\pi^{0}\to\gamma\gamma$ in chiral perturbation theory. We confirm previous estimates of this and a related LEC and obtain with additional input from theory and experiment the value $\Gamma(\pi^{0}\to\gamma\gamma)=8.07(10)\leavevmode\nobreak\ \mathrm{eV}$ in perfect agreement with recent results in the literature Kampf:2009tk . The improvements of our lattice calculation allow for a model-independent determination of the pion-pole contribution to hadronic light-by-light scattering in the muon $(g-2)$, expected to be numerically dominant, along with the $\eta$ and $\eta^{\prime}$ pseudoscalar mesons. Our final lattice result reads $a_{\mu}^{\mathrm{HLbL};\pi^{0}}=(59.7\pm 3.6)\times 10^{-11}$ which corresponds to a precision of $6\%$. With the normalization of the form factor obtained from the measurement of the decay width $\pi^{0}\to\gamma\gamma$ by the PrimEx experiment, we get the lattice and data-driven estimate $a_{\mu}^{\mathrm{HLbL};\pi^{0}}=(62.3\pm 2.3)\times 10^{-11}$ with a precision of 4%.

This paper is organized as follows. In Sec. II, we describe our methodology to extract the pion TFF. In particular, we explain how $\mathcal{O}(a)$-improvement is implemented. In Sec. III, we present our results, extrapolated to the physical point, and discuss potential sources of systematic errors. Then, in Sec. IV, after comparing our data with various phenomenological models, we use our result to compute relevant phenomenological quantities, including the pion decay $\pi^{0}\to\gamma\gamma$ and the pion-pole contribution in the hadronic light-by-light scattering contribution to the muon $(g-2)$. We conclude in Sec. V with a summary of our work and present some possible improvements for future calculations.

II Methodology

II.1 Extraction of the transition form factor

In this section, we use the same notations as in Ref. Gerardin:2016cqj and recall only the main equations. In Minkowski spacetime, the TFF describing the interaction between a neutral pion with momentum $p$ and two off-shell photons with momenta $q_{1}$ and $q_{2}$ is defined via the following matrix element,

[TABLE]

where $J_{\mu}$ is the hadronic component of the electromagnetic current, $p=q_{1}+q_{2}$ and $\epsilon_{\mu\nu\alpha\beta}$ is the fully antisymmetric tensor with111Latin indices run from 1 to 3 and Greek indices run from 0 to 3 $\epsilon^{0123}=1$. In Euclidean spacetime, the TFF is obtained after analytical continuation. The latter is valid for $q_{1,2}^{2}<s_{0}$, with $s_{0}$ the threshold for hadron production in the vector channel Ji:2001wha ; Ji:2001nf 222In the isovector case, the threshold $s_{0}$ is given by $4m_{\pi}^{2}$ or the square $\rho$ meson mass, depending on how light the quarks are., and reads

[TABLE]

where $n_{0}$ denotes the number of temporal indices carried by the two vector currents, $\omega_{1}$ is a real free parameter such that $q_{1}=(\omega_{1},\vec{q}_{1})$ and the superscript E stands for Euclidean. It is convenient to write Eq. (2) as

[TABLE]

where $E_{\pi}$ is the pion energy, $Z_{\pi}=\langle 0|P(0)|\pi\rangle$ is the overlap of the pseudoscalar operator with the pion state333Fixing the phase of the pion state via the relation $\langle 0|A_{\mu}^{a}(x)|\pi^{b},\vec{p}\rangle=iF_{\pi}p_{\mu}\delta^{ab}e^{-ip\cdot x}$, with $A_{\mu}^{a}=\bar{\psi}\gamma_{\mu}\gamma_{5}\frac{\tau^{a}}{2}\psi$ the isovector axial current, the partially conserved axial current (PCAC) relation implies that the overlap $Z_{\pi}=-\sqrt{2}im_{\pi}^{2}F_{\pi}/m$ of the operator $P=\bar{\psi}\gamma_{5}\frac{\tau^{3}}{\sqrt{2}}\psi$ is purely imaginary. Here $F_{\pi}\simeq 92\,$MeV is the pion decay constant and $m$ is the average up/down quark mass. and $\widetilde{A}_{\mu\nu}(\tau)$ is related to the three-point correlation function computed on the lattice,

[TABLE]

The three-point correlation function is defined as

[TABLE]

where $\tau=t_{i}-t_{f}$ is the time separation between the two vector currents and $t_{\pi}={\rm min}(t_{f}-t_{0},t_{i}-t_{0})$ is the minimal time separation between the pion interpolating operator and the two vector currents. Finally, even if not explicitly written, the functions $\widetilde{A}_{\mu\nu}(\tau)$, as well as the three-point correlation functions, depend on the photon spatial momenta $\vec{q}_{1}$ and $\vec{q}_{2}$.

II.2 Orbits and kinematical reach

The TFF depends on the two virtualities $q_{1}^{2}$ and $q_{2}^{2}$. Given that we use periodic boundary conditions in space, the kinematical range accessible on the lattice can be parametrized by

[TABLE]

Let the spatial momentum $\vec{p}$ of the pion be fixed. For a given $\vec{q}_{1}$, the second photon momentum $\vec{q}_{2}=\vec{p}-\vec{q}_{1}$ is fully determined and the system (6) describes a single curve in the $q_{1}^{2},q_{2}^{2}$ plane, parametrized by $\omega_{1}$. If we now consider a fixed value of $|\vec{q}_{1}|$, there is a finite number of realizations of $\vec{q}_{1}$, forming an orbit or a set of orbits444For instance, the vectors $\vec{n}=(3,0,0)$ and $\vec{n}=(2,2,1)$, which have the same norm, do not belong to the same lattice orbit. on the reciprocal cubic lattice. To this set of $\vec{q}_{1}$-values correspond in general several values of $|\vec{q}_{2}|^{2}$, each associated with a separate curve in the $(q_{1}^{2},q_{2}^{2})$ plane. In our calculation, the three-point function evaluations obtained for a given set of $(|\vec{q}_{1}|,\,|\vec{q}_{2}|)$ are averaged over in order to increase the statistical precision. The presence of hypercubic artifacts, due to the breakdown of rotational symmetry on the lattice, is discussed in Sec. III.4.1.

In Ref. Gerardin:2016cqj , we chose the pion to be at rest, and the corresponding orbits for the ensemble N200, whose parameters are given in Table 1, are shown in the left panel of Fig. 1. This setup is well suited to probe large virtualities in the double-virtual case, the kinematical region where no experimental data are available. However, because of their large eccentricity, these orbits are limited to rather low virtualities in the single-virtual case. Here we therefore include a second frame where the pion has one unit of momentum in the $z$-direction. In this case, we can access larger virtualities for the single-virtual form factor, as can be seen in the right panel of Fig. 1. In principle, even larger virtualities could be reached by increasing the pion momentum, but the signal-to-noise ratio would deteriorate rapidly. Moreover, since within our computational setup the pion interpolating operator is implemented using a sequential quark propagator, every new pion momentum $\vec{p}$ requires a new inversion of the Dirac operator, the most expensive part of the numerical simulation.

II.3 Decomposition of the integrand $\widetilde{A}_{\mu\nu}(\tau)$

The main quantity of interest in the lattice calculation is the function $\widetilde{A}_{\mu\nu}(\tau)$ defined through Eq. (4) and directly related to the three-point correlation function. Starting from Eq. (3), which holds for all real values of $\omega_{1}$ such that $q_{1,2}^{2}<s_{0}$, we consider the analytic continuation for complex values of $\omega_{1}=i\widetilde{\omega}$,

[TABLE]

Again, we do not write the dependence on the spatial momenta explicitly: $M_{\mu\nu}^{\rm E}$ is a function of $q_{1}$ and $q_{2}$ and depends implicitly on $\widetilde{\omega}$. Relation (7) is valid as long as $\widetilde{A}_{\mu\nu}(\tau)$ falls off exponentially for $\tau\to-\infty$ 555For $\tau\to+\infty$, an exponential falloff is guaranteed, irrespective of the kinematics.. If $E_{\pi}<\sqrt{s_{0}}$, which holds in all cases considered in this paper, the relation is guaranteed to hold for all values of $(\vec{q}_{1},\vec{q}_{2})$. The inversion of the Fourier transform then yields

[TABLE]

From Eq. (1), we have $M_{\mu\nu}=q_{\mu\nu}\,{\cal F}_{\pi^{0}\gamma^{*}\gamma^{*}}$, with

[TABLE]

where the coefficients $P_{\mu\nu}$ and $Q_{\mu\nu}$ do not depend on $\omega_{1}$. Thus, we can write

[TABLE]

with $P^{E}_{\mu\nu}=iP_{\mu\nu}$ and $Q^{E}_{\mu\nu}=(-i)^{n_{0}}Q_{\mu\nu}$ and where all the information is encoded into a single rotationally invariant function $\widetilde{A}^{(1)}$ defined through

[TABLE]

In Eq. (11) it is understood that the arguments of the TFF are given by Eq. (6) with $\omega_{1}$ set to $i\widetilde{\omega}$. In particular, the integral transform does not keep either virtuality fixed and $q_{2}^{2}$ even carries an imaginary part.

It is easiest to interpret the components of $\widetilde{A}_{\mu\nu}(\tau)$ in a spatially covariant notation using the Euclidean metric. With $\vec{\epsilon}$ and $\vec{\epsilon}^{\,\prime}$ two unit vectors, we can write

[TABLE]

These are the master relations expressing $\widetilde{A}_{\mu\nu}(\tau)$ in terms of an integral transform of the pion transition form factor. In the pion rest frame, $\widetilde{A}_{kl}(\tau)$ measures the transform $\widetilde{A}_{1}(\tau)$ of the transition form factor. At nonzero $\vec{p}$, $\widetilde{A}_{kl}(\tau)$ measures a linear combination of $\widetilde{A}_{1}(\tau)$ and its temporal derivative; in addition, the components $\widetilde{A}_{0k}(\tau)$ are proportional to $\widetilde{A}_{1}(\tau)$. For convenience, with a pion carrying one unit of momentum in the $z$-direction, we use the notation $\widetilde{A}^{(2)}(\tau)$ to write

[TABLE]

We remark that the matrix elements $M_{\mu\nu}$, being proportional to the TFF, are real. Therefore, $M^{E}_{kl}$ are real and $M^{E}_{0k}$ are imaginary. From Eq. (3), we conclude that $\widetilde{A}_{kl}$ are imaginary and $\widetilde{A}_{0k}$ are real666Remember that the pion overlap factor $Z_{\pi}$ is imaginary.. Thus, the two scalar functions $\widetilde{A}^{(1)}(\tau)$ and $\widetilde{A}^{(2)}(\tau)$ are real. For the vector meson dominance (VMD) and the lowest meson dominance (LMD) Moussallam:1994xp ; Knecht:1999gb models, explicit expressions for $\widetilde{A}_{\mu\nu}(\tau)$ are given in Appendix A.

Finally, from Eqs. (8) and (9), we deduce the symmetry

[TABLE]

which corresponds to the Bose symmetry of the two photons coupling to the pion. Here we have written explicitly the dependence of $\widetilde{A}_{\mu\nu}(\tau)$ on the spatial momenta. Equivalently, we have $\widetilde{A}^{(1)}(\tau;\vec{q}_{1},\vec{q}_{2})=\widetilde{A}^{(1)}(-\tau;\vec{q}_{2},\vec{q}_{1})e^{-E_{\pi}\tau}$. The relation (14) can be checked explicitly in the case of the VMD and LMD models using the expressions provided in Appendix A. These symmetries are exploited in the analysis of the numerical data.

II.4 Numerical integration

In the moving frame with $\vec{p}=(2\pi/L)\hat{z}$, both functions $\widetilde{A}^{(1)}(\tau)$ and $\widetilde{A}^{(2)}(\tau)$ are plotted in Fig. 2. The function $\widetilde{A}^{(1)}(\tau)$ is always positive with a cusp at $\tau=0$, related to the OPE at short distances, as discussed in Ref. Gerardin:2016cqj . On the contrary, the sign of $\widetilde{A}^{(2)}(\tau)$ changes at $\tau=0$ where the function is discontinuous. Therefore, the numerical integration is more challenging since large cancellations can occur in applying Eq. (3) to obtain the TFF. In particular, a naive replacement of the integral by a sum over discrete $\tau$ (trapezoidal rule) leads to noticeable numerical errors. A Simpson integration, however, reduces significantly this source of error as can be seen in Fig. 3 : here, the results correspond to fictitious data where the three-point correlation function is generated assuming an LMD model with parameters obtained from a fit to the ensemble N200. When $\vec{p}=\vec{0}$, the function $\widetilde{A}^{(1)}(\tau)$ is always positive and we do not observe any significant difference between the two integration schemes.

Due to the finite time extent of the lattice, and because the signal deteriorates at large time separations $\tau$, the integration in Eq. (3) cannot be performed up to infinity. However, the VMD model is expected to give a good description of the data in a wide $\tau$ window, as shown in Appendix A. Assuming a VMD parametrization, the functions $\widetilde{A}^{(1)}(\tau)$ and $\widetilde{A}^{(2)}(\tau)$, which depend on two parameters $\alpha$ and $M_{V}$, can be computed explicitly from Eqs. (11, 12a, 12b) and the results are given in Appendix A. As in Ref. Gerardin:2016cqj , we start by fitting our data using the VMD parametrization at large $\tau$ where the normalization of the TFF, $\alpha$, is treated as a free fit parameter. The vector meson mass $M_{V}$ is set to the value extracted from the vector two-point correlation function evaluated on the same ensembles but with higher statistics. This fit is a global fit where all momenta and both $\widetilde{A}^{(1)}(\tau)$ and $\widetilde{A}^{(2)}(\tau)$ are fitted simultaneously. Then, we integrate over lattice data up to some cutoff $\tau_{c}$ and use the fit to evaluate the remaining part of the integral with $|\tau|>\tau_{c}$. The value of $\tau_{c}\gtrsim 1.5\leavevmode\nobreak\ \mathrm{fm}$ is chosen sufficiently large such that the two-point vector correlation function is dominated by the ground state. In practice, we only include points where more than $80\%$ of the integrand comes from lattice data. We have also repeated the analysis with the LMD model and the difference between the VMD and LMD models is used to estimate the systematic uncertainty associated with the treatment of the tail. This systematic error is always smaller than the statistical error. Typical fits for our lightest ensemble, D200, are depicted in Fig. 2.

II.5 Lattice ensembles, correlation functions and $\mathcal{O}(a)$-improvement

This work is based on a subset of the $N_{f}=2+1$ Coordinated Lattice Simulations (CLS) ensembles Bruno:2014jqa generated using the openQCD suite Luscher:2012av . They use $\mathcal{O}(a)$-improved Wilson fermions, with the nonperturbative coefficient $c_{\rm SW}$ determined in Ref. Bulava:2013cta , and tree-level $\mathcal{O}(a^{2})$-improved Lüscher-Weisz action for the gauge field. To avoid the freezing of the topological charge when approaching the continuum limit Luscher:2012av ; Luscher:2011kk , CLS ensembles use periodic boundary conditions in space and open boundary conditions in time. The parameters of the simulations are summarized in Table 1. We use four different lattice spacings in the range 0.050 - 0.086 fm and several pion masses down to 200 MeV to perform the continuum and chiral extrapolations. For the latter, the ensembles considered here have been generated keeping the average bare quark mass $m_{\rm q}^{\rm av}=(2m_{l}+m_{s})/3$ fixed and using two degenerate light quarks $m_{u}=m_{d}=m_{l}$. This chiral trajectory has the advantage of automatically keeping the $\mathcal{O}(a)$-improved bare coupling $\tilde{g}_{0}$ constant Bietenholz:2010jr . All our ensembles included in the final analysis satisfy $m_{\pi}L>4$ such that finite-size effects are expected to be small. The scale setting was performed in Ref. Bruno:2016plf with a precision of $1\%$ using a linear combination of the pion and kaon decay constants. Statistical errors are estimated using the jackknife procedure and we propagate errors associated with the renormalization constant of the vector current, the scale setting, the pion mass $m_{\pi}$ and decay constant $f_{\pi}$ in the chiral extrapolations.

To reduce discretization effects and obtain a shorter continuum extrapolation than in Gerardin:2016cqj , on-shell $\mathcal{O}(a)$-improvement has been implemented. In addition to the improvement of the action, it requires the implementation of the renormalized $\mathcal{O}(a)$-improved vector current. Two discretizations are used, the local $(l)$ and the point-split $(c)$ lattice vector currents,

[TABLE]

where $\lambda^{a}$ are the eight Gell-Mann matrices. With the tensor current defined as $\Sigma^{a}_{\mu\nu}(x)=-\frac{1}{2}\,\overline{\psi}(x)[\gamma_{\mu},\gamma_{\nu}]\frac{\lambda^{a}}{2}\psi(x)$, the improved vector current is given by

[TABLE]

where the coefficient $c_{\mathrm{V}}$ differs for the local and point-split vector currents and has been determined nonperturbatively in Ref. Gerardin:2018kpy . The point-split vector current is conserved on the lattice and does not need to be further renormalized. For the local vector current, the renormalization pattern reads Bhattacharya:2005rb

[TABLE]

where $Z_{\mathrm{V}}$, $b_{\rm V}$ and $\overline{b}_{\rm V}$ have been evaluated nonperturbatively in Ref. Gerardin:2018kpy and $M_{\rm q}=\mathrm{diag}(m_{l},m_{l},m_{s})$ is the bare (subtracted) quark mass matrix. The renormalized coupling is given by $\widetilde{g}_{0}^{2}=g_{0}^{2}(1+b_{\rm g}am_{\rm q}^{\rm av})$. The coefficient $f_{\rm V}$ starts at order $\mathcal{O}(g^{6})$ in perturbation theory and is neglected here. For the electromagnetic current $J_{\mu}$ with up, down and strange quarks, it is convenient to use the isospin decomposition

[TABLE]

where $V_{\mu}^{a}=\bar{\psi}\gamma_{\mu}\frac{\lambda^{a}}{2}\psi$ is the octet of vector currents and the hat means that the current is both improved and renormalized. In particular, Eq. (17) reduces to

[TABLE]

with $V^{0,I}_{\mu}=\frac{1}{2}\bar{\psi}\gamma_{\mu}\psi$ the flavor-singlet current and

[TABLE]

We use $P(x)=(1/\sqrt{2})(\overline{u}(x)\gamma_{5}u(x)-\overline{d}(x)\gamma_{5}d(x))$ as an interpolating operator for the neutral pion. The renormalized and $\mathcal{O}(a)$-improved three-point correlation function appearing in Eq. (5), obtained with two local vector currents, reads

[TABLE]

where the first line corresponds to the connected part and the second and third lines to the disconnected parts. We did not explicitly write the improvement term proportional to $c_{\mathrm{V}}$. Here, $G_{l}$ and $G_{s}$ denote the light and strange quark propagators respectively. If $f_{\rm V}$ is neglected, the connected part renormalizes proportionally to $Z_{3}^{2}$, since $Z_{8}+\sqrt{3}Z_{80}=Z_{3}-2Z_{\mathrm{V}}f_{\rm V}\,a(m_{s}-m_{l})$. Similarly, with one local and one conserved vector currents, one finds

[TABLE]

where a summation over $y^{\prime}$ and $y^{\prime\prime}$ is understood and where

[TABLE]

The lattice quark propagator is $\gamma_{5}$-Hermitian, $G^{{\dagger}}=\gamma_{5}G\gamma_{5}$, and the linear operator $\mathcal{V}^{(y)}_{\mu}$ is anti $\gamma_{5}$-Hermitian, $\mathcal{V}^{(y){\dagger}}_{\mu}=-\gamma_{5}\mathcal{V}^{(y)}_{\mu}\gamma_{5}$.

The three-point correlation function is computed using the same technique as in Ref. Gerardin:2016cqj . A point source is created on the time slice $y_{0}$ with a random value of $\vec{y}$ and a sequential propagator is computed at $x_{0}$, for two values of the pion spatial momenta $\vec{p}=\vec{0}$ and $\vec{p}=(2\pi/L)\hat{z}$. With our computational setup we obtain all values of the photon momenta $\vec{q}_{1}$ and $\vec{q}_{2}$ such that $\vec{q}_{1}+\vec{q}_{2}=\vec{p}$ without any new inversion of the Dirac operator. The number of sources per gauge configuration is given in Table 1 and the location of $y_{0}$ is discussed in Sec. III.1. For all ensembles, we are able to probe photon virtualities up to $3\leavevmode\nobreak\ \mathrm{GeV}^{2}$ in the double-virtual case and up to $1.5\leavevmode\nobreak\ \mathrm{GeV}^{2}$ in the single-virtual case. This is a major improvement compared to our previous study, where we were limited to virtualities below $0.5\leavevmode\nobreak\ \mathrm{GeV}^{2}$ in the single-virtual case.

In Eq. (16), we have some freedom in the choice of the lattice derivative: different definitions differ by $\mathcal{O}(a)$ effects and therefore contribute only to order $a^{2}$ in correlation functions. We have considered the symmetric, the backward and the forward derivatives, respectively, given by

[TABLE]

For the vector current located at the source, we choose the forward derivative $\partial_{\mu}^{f}f(x)$, since a symmetric derivative would require more inversions of the Dirac operator. For the vector current located at the sink, we have compared the three different discretizations. We define the ratio

[TABLE]

where $\widetilde{A}^{(1),{\rm lc}}(\tau)$ and $\widetilde{A}^{(1),{\rm ll}}(\tau)$, defined in Eq. (3), are, respectively, computed using one or two local vector currents and $\delta\widetilde{A}^{(1),{\rm lc}}(\tau)$ is the statistical error associated with $\widetilde{A}^{(1),{\rm lc}}(\tau)$. In the left panel of Fig. 4, this ratio is plotted with and without improvement, using different discretizations of the lattice derivative. After $\mathcal{O}(a)$-improvement, the discrepancy between both discretizations is reduced as expected and the symmetric derivative turns out to be the best choice.

Finally, we note that on-shell improvement of correlation functions is not sufficient for $\tau={\rm O}(a)$. For $\tau=\pm a$, we do not use the symmetric derivative but the forward and backward derivatives, respectively, to avoid using the correlation function at $\tau=0$ where $A^{(2)}(\tau)$ is discontinuous (see right panel of Fig. 2)). In the right panel of Fig. 4 we plot $\widetilde{A}(\tau)$ using both discretizations ($\widetilde{A}^{ll}(\tau)$ and $\widetilde{A}^{lc}(\tau)$) with and without improvement of the vector currents using the symmetric derivative. We observe that after improvement, with the exception of $\tau=0$, the two discretizations are in excellent agreement with each other.

III Results

III.1 Two-point pseudoscalar correlation function

The pion energy $E_{\pi}(\vec{p})$ and its overlap $Z_{\pi}$ with our interpolating operator are extracted from the pseudoscalar two-point correlation function. In the limit of large source-sink time separation, and with both source and sink far from the boundary, the correlation function behaves as

[TABLE]

We want to extract $E_{\pi}(\vec{p})$ and $Z_{\pi}$ in a region where the excited-state contribution and boundary effects are both small. We therefore place the source far from the boundary located at $x_{0}=0$. In Fig. 5, we show the result for the effective mass defined by

[TABLE]

For large source-sink time separation, we observe a plateau up to some time where the sink starts to be close to the boundary located at $x_{0}=T-a$. In the case of open boundary conditions in time, and including only the first excited state, we first fit our data using the Ansatz

[TABLE]

where $E_{\pi}^{\prime}(\vec{p})$ and $Z_{\pi}^{\prime}$ are the energy and the overlap of the first excited state with momentum $\vec{p}$ in the pseudoscalar channel and the third term includes the first boundary excited state: a finite volume two-pion state with vanishing momentum. In the fits, we assume that $E_{2\pi}=2E_{\pi}(\vec{0})$. Explicitly, denoting $|B\rangle$ the boundary state, the prefactor $A(\vec{p})$ is given by

[TABLE]

To reduce the number of fit parameters, both momenta are fitted simultaneously. In a second step, the data are fitted using a single-exponential ansatz in the region where excited state contributions are small compared to the statistical precision, and the results are given in Table 2. We note that the pion energies are compatible with the relativistic dispersion relation. Finally, we have checked that the source is sufficiently far from the boundary such that the overlaps are correctly extracted: the last term in Eq. (28) contributes to $Z_{\pi}$ as less than half its statistical error. The same value of $y_{0}$ is used to compute the three-point correlation function. The statistical error on the ratio $E_{\pi}/Z_{\pi}$ in Eq. (3) lies between 0.5% and 0.8%.

The chiral extrapolation of the TFF is done using the dimensionless parameter $\widetilde{y}=m_{\pi}^{2}/(16\pi^{2}f_{\pi}^{2})$. In order to propagate the error on $m_{\pi}$ and $f_{\pi}$, we also compute the pion decay constant on the same gauge configurations for each ensemble. The bare matrix element of interest is extracted from

[TABLE]

where the axial-pseudoscalar two-point correlation function is given by

[TABLE]

The $\mathcal{O}(a)$-improved axial current reads $A^{I}_{0}(x)=A_{0}(x)+c_{\rm A}\,\partial^{(s)}_{0}P(x)$ where the improvement parameter $c_{\rm A}$ has been determined nonperturbatively in Ref. Bulava:2015bxa . The ratio $R(x_{0},y_{0})$ is fitted using the Ansatz

[TABLE]

From this bare matrix element, the renormalized and $\mathcal{O}(a)$-improved pseudoscalar decay constant is given by777Within our conventions, the physical value of $f_{\pi}$ is $F_{\pi}=92.4\leavevmode\nobreak\ \mathrm{MeV}$. We distinguish between $F$ the pion decay constant in the chiral limit, $F_{\pi}$ the pion decay constant at physical pion mass and $f_{\pi}$ the pion decay constant at unphysical pion mass.

[TABLE]

where $Z_{\rm A}$ is the renormalization factor of the axial current and $b_{\rm A}$, $\overline{b}_{\rm A}$ are improvement parameters DallaBrida:2018tpn ; Korcyl:2016ugy . The results are summarized in Table 2. In practice, we correct these results for finite-size effects (FSE) using chiral perturbation theory ($\chi$PT) as described in Ref. Colangelo:2005gd . These corrections are small and we find that they correctly account for FSE on the ensembles H105/N101, which were generated using the same action parameters but different lattice volumes. However the FSE correction fails on ensemble H200, which corresponds to our smallest volume ($L\approx 2$ fm) and was not used to extract the pion TFF.

As a consistency check, we have performed an extrapolation to the physical point assuming the functional form

[TABLE]

inspired from NLO $\chi PT$ Gasser:1983yg and where we allow for a quadratic dependence in $\tilde{y}$. The four fit parameters are $f_{\pi}(0,\tilde{y}_{\mathrm{phys}})$, $\bar{\ell}_{4}$, $\epsilon_{f_{\pi}}$ and $\delta_{f_{\pi}}$. Our result at the physical point $f_{\pi}=92.1(1.8)\leavevmode\nobreak\ \mathrm{MeV}$ and the value $\bar{\ell}_{4}=3.56(44)_{\mathrm{stat}}(80)_{\mathrm{syst}}$ turns out to be in good agreement with the Flavour Lattice Averaging Group (FLAG) average Aoki:2019cca .

III.2 The transition form factor

The results for the TFF at our lightest pion mass, corresponding to ensemble D200, are plotted in Fig. 6 for two different kinematics : the single-virtual case and the double-virtual case with $Q_{1}^{2}=Q_{2}^{2}$. Black and blue points correspond to the TFF obtained in the pion rest frame and moving frame respectively.

In the pion rest frame, results obtained using $A^{(1)}(\tau)$ are statistically more precise when both photons carry the same virtuality (right panel of Fig. 6). In this case $\omega_{1}=m_{\pi}/2$ and the full integrand in Eq. (3) is symmetric. In the single-virtual case, $\omega_{1}$ takes larger values and the integral probes further the tail of the function $\widetilde{A}^{(1)}(\tau)$, where the noise over signal ratio increases rapidly with $\tau$: the signal is lost for virtualities above $0.5\leavevmode\nobreak\ \mathrm{GeV}^{2}$. In the moving frame, the situation is similar for $\widetilde{A}^{(1)}(\tau)$ (even if the integrand is not exactly symmetric in the double-virtual case, unless $|\vec{q}_{1}|^{2}\approx|\vec{q}_{2}|^{2}$). For $\widetilde{A}^{(2)}(\tau)$, since the sign of the function changes at $\tau=0$, the situation is the opposite: results are more precise in the single-virtual case where there are fewer cancellations between positive and negative contributions and we can reach higher virtualities (see right panel of Fig. 2).

We also point out that results obtained in both frames are fully consistent, confirming that sources of breaking the Lorentz invariance in the lattice calculation are small. In particular, since the two frames are affected by different sources of finite volume effects, it is a first hint that FSE are small. A more detailed study of FSE is done in Sec. III.4.2.

III.3 Parametrization of ${\cal F}_{\pi^{0}\gamma^{*}\gamma^{*}}$ and extrapolation to the physical point

In this section we propose a method to extract the TFF over the whole kinematical range, and at the physical point, based on its analytical properties. We introduce the conformal variables $z_{1}$ and $z_{2}$ defined through Boyd:1995sq

[TABLE]

which map the branch cut, starting at $t_{c}=4m_{\pi}^{2}$, onto the unit circle $|z_{k}|=1$. Here, $t_{0}$ is a free parameter representing the virtuality mapping to $z_{k}=0$. Since the TFF is analytic for $|z_{k}|<1$, one can write

[TABLE]

where the coefficients $c_{nm}=c_{mn}$ are symmetric due to the Bose symmetry. Since $|z_{k}|<1$ one can expect a fast convergence of the sum, where only a few terms with $n,m\leq N$ are needed at a given accuracy. The optimal choice for $t_{0}$, which reduces the maximum value of $|z_{k}|$ in the range $[0,Q_{\max}^{2}]$, is given by

[TABLE]

This maximum is reached at $z_{k}=0$ and $z_{k}=Q^{2}_{\max}$. Using the values $m_{\pi}=134.9\leavevmode\nobreak\ \mathrm{MeV}$ and $Q_{\max}^{2}=4\leavevmode\nobreak\ \mathrm{GeV}^{2}$, one finds $|z_{\max}|=0.46$ well below one. The coefficients $c_{nm}$ satisfy the relation

[TABLE]

which ensures that the coefficients $c_{nm}$ are not only bounded but also decrease to zero for sufficiently large $n,m$. In practice, the TFF can be multiplied by any analytical function $P(Q_{1}^{2},Q_{2}^{2})$ and the resulting product expanded in powers of the $z_{k}$. Some choices may improve the convergence rate of the series expansion. At short distances, the behavior of the TFF is predicted by the Brodsky-Lepage behavior in the single-virtual case and by the OPE in the double-virtual case BL_3_papers ; Nesterenko:1982dn ; Novikov:1983jt ,

[TABLE]

Therefore, it is convenient to consider

[TABLE]

where $M_{V}=775\leavevmode\nobreak\ \mathrm{MeV}$ is the vector meson mass. With this choice, the parametrization of the TFF at the finite value of $N$ decreases asymptotically as $1/Q^{2}$ in all directions in the $(Q_{1}^{2},Q_{2}^{2})$ plane, in accord with the Brodsky-Lepage and the OPE behavior. This feature is one major reason for us preferring the $z$-expansion over the LMD+V model Knecht:2001xc used in Gerardin:2016cqj , which e.g. for constant $Q_{1}^{2}\neq 0$ and $Q_{2}^{2}\to\infty$ does not vanish, even though the corresponding integrals for the pion-pole contribution to the hadronic light-by-light scattering in the muon $g-2$ are still convergent with the LMD+V model (see Sec. IV.4). However, at the precision level we aim for, the proper high-energy behavior in all directions is important. It should be noted that we do not impose the explicit values of the coefficients on the right-hand side of Eqs. (39) and (40), since they receive higher-order corrections in perturbative QCD and at higher twist, see the discussion below.

The imaginary part of the TFF behaves as $(q^{2}-t_{c})^{3/2}$ near threshold (P-wave). This property is not fulfilled at finite $N$ but888The authors thank Martin Hoferichter for pointing out this fact to us., as shown in Ref. Bourrely:2008za , this constraint can be implemented by imposing

[TABLE]

which leads to the modified $z$-expansion

[TABLE]

Finally, to take into account the discretization effects and the dependence on the quark masses used in the simulations, the coefficients $c_{nm}$ are allowed to vary linearly with the variable $\widetilde{y}=m_{\pi}^{2}/(16\pi^{2}f_{\pi}^{2})$ and quadratically with the lattice spacing,

[TABLE]

where $d=1,2$ stands for the two discretizations of the three-point correlation function. We perform a global fit and the results are summarized in Table 3. Due to the large set of data, the correlation matrices are ill-conditioned and we perform uncorrelated fits. The error on the fit parameters is then obtained from a jackknife procedure, using blocking to take into account autocorrelations. The results for $N=3$ are shown in Fig. 7. The chi-square by degree of freedom is $\chi^{2}/\mathrm{d.o.f.}=1.1$ and the correlation matrix of the coefficients in the continuum and at the physical pion mass is

[TABLE]

In Appendix B, we study the systematic error associated with the truncation of the sum in Eq. (43) and we conclude that $N=3$ is already sufficient for an accuracy at the percent level. We also checked that using $P(Q_{1}^{2},Q_{2}^{2})=1$ leads to compatible results in the range where we have lattice data. In the single-virtual case, our results are in good agreement with experimental data. In the single-virtual case and in the double-virtual case, we also find good agreement with the results obtained in the dispersive framework Hoferichter:2018dmo ; see Fig. 7. A more detailed comparison with phenomenology is provided in the next section.

Finally, we remark that a different fit strategy would consist in fitting each ensemble independently and then extrapolating each coefficient $c_{nm}$ to the physical point. First, we found that the global fit procedure is more stable. Second, the continuum and chiral extrapolation of the individual coefficients $c_{nm}$ is nontrivial, as they are correlated. Nonetheless, we provide the $z$-expansion on three individual lattice ensembles in Appendix C to allow for comparisons prior to the chiral and continuum extrapolation.

III.4 Systematic errors

III.4.1 Hypercubic effects

At finite lattice spacing, the O(3) rotational symmetry is broken down to the cubic subgroup H(3). Spatial momenta, equivalent in the continuum, but belonging to different H(3) orbits may be affected by different lattice artifacts. This has been observed in previous calculations of form factors (recently in Ref. Lubicz:2017syv for example). As explained in Sec. II.2, in the calculation of the TFF, we have averaged over all equivalent combinations of $(\vec{q}_{1}^{\ 2},\vec{q}_{2}^{\ 2})$. In this section, we study the validity of this approach.

Any polynomial function of the spatial momenta $\vec{q}$, and invariant under H(3), can be expressed in terms of the three invariants

[TABLE]

Therefore, each H(3) orbit is uniquely characterized by the values of $q_{i}^{[2]}$, $q_{i}^{[4]}$ and $q_{i}^{[6]}$ with $i=1,2$.

In the pion rest frame, both photons have back-to-back spatial momenta and $\vec{q}_{1}=-\vec{q}_{2}$. The first kinematical configuration with two orbits corresponds to $|\vec{q}_{1}|^{2}=9(\frac{2\pi}{L})^{2}$ with, for example, $\vec{q}_{1}=\frac{2\pi}{L}(\pm 3,0,0)$ (6 elements) and $\vec{q}_{1}=\frac{2\pi}{L}(\pm 2,\pm 2,\pm 1)$ (24 elements). When the pion has one unit of momentum in the $\hat{z}$-direction, equivalent $(q_{1}^{2},q_{2}^{2})$ can also be obtained from different values of $\vec{q}_{1}$ not related by H(3) symmetry but by the permutation of $q_{1}$ and $q_{2}$.

In Fig. 8, we show the results for the transition form factor without averaging over equivalent momenta. We observe that, at our level of statistical precision, we are not yet sensitive to hypercubic artifacts.

III.4.2 Finite-size effects

Two sets of ensembles (H200/N202 and H105/N101, listed in Table 1) have been generated using the same bare lattice parameters $\kappa$ and $\beta$ but with different volumes ($L/a=32$ and $48$). The ensembles N202 and N101 correspond to large volumes with $m_{\pi}L\geq 6.5$ where FSE are expected to be negligible. The two other ensembles satisfy $m_{\pi}L\approx 4.0$. The pion decay constant has been computed on all ensembles (see Table 2) with a statistical precision below $1.0\leavevmode\nobreak\ \%$ and finite-size effects are sizable: we observe a 2-3 $\sigma$ discrepancy between H200/N202 and FSE corrections using $\chi$PT Colangelo:2005gd failed to explain the difference. However, for H105 and N101, the results are in perfect agreement after FSE corrections. Since finite-size effects strongly depend on the observable (and on the estimator), the TFF has been computed on these two sets of ensembles. At our level of statistical precision, we do not observe any significant differences between different volumes and the results for H200, compared to N202, are shown in Fig. 9. Since H200 is our smallest ensemble, we neglect finite-size effects in the following and exclude the ensembles H200 from the analysis. To be conservative, we also exclude H105.

III.4.3 Quark-disconnected contributions

The quark-disconnected contributions involve a single quark loop which couples to the vector current. Within the Mainz group, they have been coded MainzDisc and computed as part of various physics projects Ce:2018ziv ; Djukanovic:2018iir ; Gerardin:2018sin on a number of ensembles (H105, N203, N200, D200 and N302). These contributions vanish exactly for the three ensembles at the heaviest pion mass, where $m_{l}=m_{s}$. The single-propagator loops have been computed using two random sources of 512 hierarchical probing vectors Stathopoulos:2013aci per configuration and for all spatial momenta with $|\vec{q}|^{2}\leq 12\times(\frac{2\pi}{L})^{2}$. The two-point correlation functions between a pseudoscalar and a vector current have been computed using stochastic sources with U(1) noise as in Ref. Gerardin:2016cqj . We define the ratio

[TABLE]

where ${\cal F}_{\pi^{0}\gamma^{*}\gamma^{*}}^{\rm conn}+{\cal F}_{\pi^{0}\gamma^{*}\gamma^{*}}^{\rm disc}$ and ${\cal F}_{\pi^{0}\gamma^{*}\gamma^{*}}^{\rm conn}$ correspond to the TFF obtained with or without including the disconnected contribution. The results are shown in the left panel of Fig. 10. We remark that the tail of the disconnected correlator was treated in the same manner as in the connected correlator. As expected, the disconnected contribution increases as we approach the physical pion mass: the CLS ensembles used in this work were generated using a constant value of the trace of the bare quark masses. Comparing ensembles generated at the same pion mass but with different lattice spacings, we do not observe significant discretization effects. For our ensemble with the lightest pion (D200), the TFF with and without including the disconnected contribution is depicted in the right panel of Fig. 10 : at our level of statistical precision, the disconnected contribution can be neglected. The impact of the disconnected contribution on the pion-pole contribution to the muon $g-2$ is discussed in more detail in Sec IV.4.3.

IV Phenomenological applications

IV.1 Parametrizing and comparing the TFF with phenomenological models

As in Ref. Gerardin:2016cqj , we compare our results with phenomenological models which have been applied to the pion-pole contribution to hadronic light-by-light scattering in the muon $g-2$. In particular, we consider the VMD model, the LMD model Moussallam:1994xp ; Knecht:1999gb and the LMD+V model Knecht:2001xc . These models are parametrized by

[TABLE]

and we refer the reader to Ref. Gerardin:2016cqj for a more detailed description of the models and phenomenological values of the parameters. We just recall the main properties relevant for this study. All models satisfy the anomaly constraint at vanishing momenta Adler:1969gk ; Bell:1969ts

[TABLE]

by setting $\alpha=\alpha^{\rm th}=1/(4\pi^{2}F_{\pi})=0.274\leavevmode\nobreak\ \mathrm{GeV}^{-1}$. The VMD model does not fulfill the OPE constraint given by Eq. (40) but rather behaves as ${\cal F}_{\pi^{0}\gamma^{*}\gamma^{*}}^{\mathrm{VMD}}(-Q^{2},-Q^{2})\sim 1/Q^{4}$ at large virtualities. The LMD model has the advantage to fulfill the OPE constraint but does not satisfy the Brodsky-Lepage (BL) behavior given by Eq. (39) since ${\cal F}_{\pi^{0}\gamma^{*}\gamma^{*}}^{\mathrm{LMD}}(-Q^{2},0)\sim\beta/M_{V}^{2}$ at large virtualities. Finally, the LMD+V model satisfies the short-distance constraints both for $Q_{1}^{2}=Q_{2}^{2}\to\infty$ and for $Q_{1}^{2}=0$, $Q_{2}^{2}\to\infty$ if one sets $\widetilde{h}_{1}=0$.

We have performed global fits using the VMD, LMD and LMD+V Ansätze where the fit parameters are, respectively, $(\alpha,M_{V})$, $(\alpha,\beta,M_{V})$ and $(\alpha,\widetilde{h}_{0},\widetilde{h}_{2},\widetilde{h}_{5},M_{V_{1}},M_{V_{2}})$. As for the $z$-expansion, each fit parameter is assumed to depend linearly on the dimensionless variable $\widetilde{y}=m_{\pi}^{2}/(16\pi^{2}f_{\pi}^{2})$ and quadratically on the lattice spacing $a$. The fits are uncorrelated. Since we have computed the TFF using two different discretizations of the vector current, we actually perform a combined fit such that any parameter has the functional form

[TABLE]

In Ref. Gerardin:2016cqj , we have shown that the VMD model fails to reproduce the lattice data at large $Q^{2}$ due to the wrong asymptotic behavior in the double-virtual case. Fitting our new data, we obtain $\chi^{2}/\mathrm{d.o.f.}=4.8$. However, we obtain a reasonable $\chi^{2}$ if we restrict the fit to the single-virtual TFF and the result at the physical point reads

[TABLE]

with $\chi^{2}/\mathrm{d.o.f.}=1.2$. The value of $\alpha^{\mathrm{VMD}}$ is compatible with the theoretical prediction for the anomaly within 2$\sigma$ and $M_{V}$ is close to the rho meson mass. For the LMD model, we fit $\alpha$, $\beta$, and $M_{V}$ on the full kinematical range. The result of the global fit at the physical point reads

[TABLE]

with $\chi^{2}/\mathrm{d.o.f.}=1.5$. The value of $\alpha^{\mathrm{LMD}}$ is compatible with the anomaly constraint with a statistical precision of $2\leavevmode\nobreak\ \%$. However, the value of $\beta$ is lower, in absolute value, than the OPE prediction at leading order $\beta^{\rm OPE}=-F_{\pi}/3=-30.8\leavevmode\nobreak\ \mathrm{MeV}$. This point will be discussed in Sec. IV.3: the OPE prediction neglects higher twists and $\mathcal{O}(\alpha_{s})$ corrections which are sizable at virtualities accessible on the lattice. In Fig. 11 we show the quality of the fit for the ensembles N203 and N200 in the single-virtual configuration. At large $Q^{2}$, we observe a significant deviation between the fit and the lattice data which explains the relatively bad $\chi^{2}$. This feature was not observed in our previous work Gerardin:2016cqj where only virtualities $Q^{2}<0.5\leavevmode\nobreak\ \mathrm{GeV}^{2}$ were accessible in the single-virtual case. The fact that the LMD model fails to describe the lattice data in the single-virtual case is not unexpected since the model has the wrong asymptotic behavior. This explains the discrepancy between the TFF extrapolated to the physical point using the LMD model and the experimental data, as shown in Fig. 12.

Finally, for the LMD+V model we explicitly set $\widetilde{h}_{1}=0$ to satisfy the OPE constraint. Furthermore, to reduce the number of fit parameters, and inspired by quark models, we assume a constant shift in the spectrum and set $M_{V_{2}}(\widetilde{y})=m^{\exp}_{\rho^{\prime}}+M_{V_{1}}(\widetilde{y})-m^{\exp}_{\rho}$ with $m^{\exp}_{\rho^{\prime}}=1.465\leavevmode\nobreak\ \mathrm{GeV}$. All the remaining parameters are fitted and the result reads

[TABLE]

with $\chi^{2}/\mathrm{d.o.f.}=1.2$. We note that the same fit would be unstable if only lattice data obtained in the pion rest frame were included. The coefficient $\widetilde{h}_{0}$, related to the OPE behavior, is in good agreement with the leading order OPE prediction $\widetilde{h}_{0}=-F_{\pi}/3$. The phenomenological value $\widetilde{h}_{2}=0.327\leavevmode\nobreak\ \mathrm{GeV}^{2}$ Melnikov:2003xd can be fixed by comparing with the higher-twist corrections in the OPE in Eq. (40) and our result turns out to be in good agreement as well. Finally, $\widetilde{h}_{5}$ is consistent with the phenomenological prediction $\widetilde{h}_{5}=-0.166\pm 0.006\leavevmode\nobreak\ \mathrm{GeV}$ obtained in Ref. Knecht:2001xc by fitting the LMD+V model to the CLEO data in the single-virtual case. At our level of precision, this model provides a good description of our lattice data in the whole kinematical range covered by this study. It is also perfectly compatible with the available experimental data, as can be seen on the left panel of Fig. 12.

Recently, a new parametrization of the pion TFF over the whole spacelike region has been proposed in Danilkin:2019mhd . The parametrization of the $(Q_{1}^{2},Q_{2}^{2})$ dependence relies on a single parameter, which can be obtained by adjusting the model to the experimental data in the single-virtual configuration. The model obeys by construction the leading prediction for the asymptotic behavior at large $Q^{2}$; however, it does not accurately reproduce our lattice result at intermediate virtualities: the prediction underestimates our result by about 20% at $Q_{1}^{2}=Q_{2}^{2}=0.5\leavevmode\nobreak\ \mathrm{GeV}^{2}$ as can be seen in Fig 7. Nonetheless, in the future this model could perhaps be used as the factor $P(Q_{1}^{2},Q_{2}^{2})$ in the expansion (43).

IV.2 Normalization of the TFF, the decay width

$\Gamma(\pi^{0}\to\gamma\gamma)$ and the slope of the form factor

To leading order in the fine-structure constant $\alpha_{\rm e}$, the decay width of $\pi^{0}\to\gamma\gamma$ is given by the pion TFF as follows: $\Gamma(\pi^{0}\to\gamma\gamma)=(\pi\alpha_{\rm e}^{2}m_{\pi}^{3}/4){\cal F}_{\pi^{0}\gamma^{*}\gamma^{*}}(0,0)^{2}$. From the measurement $\Gamma(\pi^{0}\to\gamma\gamma)=7.82(22)\leavevmode\nobreak\ {\rm eV}$ by the PrimEx experiment Larin:2010kq we get the TFF at the origin with a precision of 1.4%, $\alpha^{\rm PrimEx}=0.276(4)\leavevmode\nobreak\ \mathrm{GeV}^{-1}$, and this is therefore an important benchmark for our lattice calculation. The PrimEx-II experiment aims at reducing the error on the decay width by a factor of two and a preliminary result $\Gamma(\pi^{0}\to\gamma\gamma)=7.80(13)\leavevmode\nobreak\ {\rm eV}$ has been released Gasparian:2016oyl ; Primex2Mainz , which leads to $\alpha^{\rm PrimEx-II}=0.275(2)\leavevmode\nobreak\ \mathrm{GeV}^{-1}$ with 0.8% precision. A tension exists between the state-of-the-art chiral predictions of the decay width pion_decay_NLO ; Ananthanarayan:2002kj ; pion_decay_updates ; Kampf:2009tk and the new, more precise, preliminary measurement.

In the chiral limit, the normalization of the TFF is exactly predicted by the ABJ chiral anomaly Adler:1969gk ; Bell:1969ts ,

[TABLE]

where $F$ is the pion decay constant in the chiral limit. At finite quark mass, this result receives corrections which can be computed in the framework of chiral perturbation theory ($\chi$PT) or directly using lattice QCD. The chiral expansion of the pion TFF at vanishing momenta is known up to next-to-next-to-leading (NNLO) order in $\chi$PT Kampf:2009tk . Chiral logarithms are absent at NLO when the result is expressed in terms of the physical value of $F_{\pi}$ pion_decay_NLO . They appear only at NNLO but were shown to contribute at the permille level and are negligible at our level of accuracy Kampf:2009tk . This motivates the extrapolation of the normalization of the form factor on the lattice using the Ansatz

[TABLE]

already used in a previous lattice calculation Feng:2012ck . This functional form differs from the one assumed in the previous sections but both parametrizations lead to compatible results at our level of precision. Equation. (55) has the advantage to offer the possibility to extract the low-energy constant (LEC) $C_{7}^{\mathrm{Wr}}$ in the odd-intrinsic-parity sector of chiral perturbation theory at order $p^{6}$ Bijnens:2001bb via the relation Ananthanarayan:2002kj ; Kampf:2009tk

[TABLE]

by varying the quark masses and thus $m_{\pi}$ over a range of values on our lattice ensembles. As explained in Sec. II.5, we note that our lattice ensembles lie in the chiral trajectory where the trace of the bare quark matrix is kept constant. If we repeat the analysis done in Sec. III.3 but using $f_{\pi}\,{\cal F}_{\pi^{0}\gamma^{*}\gamma^{*}}(Q_{1}^{2},Q_{2}^{2})$ instead of ${\cal F}_{\pi^{0}\gamma^{*}\gamma^{*}}(Q_{1}^{2},Q_{2}^{2})$ as our primary observable, and limiting the fit range below 1 GeV, we obtain the result

[TABLE]

where the first error is statistical and the second error is an estimate of the systematic error due to the truncation on the $z$-expansion (see Appendix B). The choice $N=1$ is already sufficient to get a $\chi/\mathrm{d.o.f.}=1.1$. The normalization of the TFF is almost unchanged compared to the fit done in Sec. III.3. The result in Table 3, with $N=3$, leads to $\alpha=0.261(13)(2)\leavevmode\nobreak\ \mathrm{GeV}^{-1}$. The value of the normalization $\alpha$ is compatible with the theoretical prediction $\alpha^{\rm th}=0.274\leavevmode\nobreak\ \mathrm{GeV}^{-1}$, given by Eq. (54) with $F$ replaced by $F_{\pi}$, with a precision of $3.5\leavevmode\nobreak\ \%$. At this level of precision, we are not yet sensitive to the other chiral corrections, not taken into account by this replacement. We note that our choice of $t_{0}$ in the $z$-expansion is optimal for the pion-pole contribution to $a_{\mu}^{\rm HLbL}$, discussed in the next subsection, but not for the normalization, as explained in Appendix B. The results for $\alpha$ and $C_{7}^{\mathrm{Wr}}$ obtained from a fit based on the LMD+V model are compatible within error bars with the model-independent results in Eq. (57).

The central value of our result for the LEC $C_{7}^{\mathrm{Wr}}$ in Eq. (57) is larger than the bound $|C_{7}^{\mathrm{Wr}}|<0.06\times 10^{-3}\leavevmode\nobreak\ \mathrm{GeV}^{-2}$ used in Ref. Kampf:2009tk . Of course with our current uncertainty we cannot exclude that this LEC vanishes. The bound is essentially based on estimates using a resonance Lagrangian with heavy pseudoscalars and the nonobservation of the decay $\pi(1300)\to\gamma\gamma$ Moussallam:1994xp ; Ananthanarayan:2002kj ; pdg:2018 . The LEC $C_{7}^{\mathrm{Wr}}$ also vanishes exactly Knecht:2001xc in simple resonance Lagrangians with vector mesons only, such as the VMD model for the TFF or the model proposed in Ref. Prades:1993ys , that do not obey short-distance constraints from the OPE. On the other hand, using the LMD model Moussallam:1994xp ; Knecht:1999gb or a resonance Lagrangian with vector and heavy pseudoscalar mesons (LMD+P) that obeys these short-distance constraints, one obtains the estimate $C_{7}^{\mathrm{Wr}}=0.35(7)\times 10^{-3}\leavevmode\nobreak\ \mathrm{GeV}^{-2}$ Kampf:2011ty . The result is actually dominated by the LMD part and the additional contribution from the heavy pseudoscalar is again estimated to be very small. Note, however, that the error might be underestimated. All resonance estimates of LECs are based on the assumption of narrow resonances in large-$N_{c}$ QCD and carry an intrinsic uncertainty of 20%-30%. Under the assumption that the LECs have a natural size of $0.5\times 10^{-3}\leavevmode\nobreak\ \mathrm{GeV}^{-2}$, one could argue that this implies rather an absolute error of $\pm(0.10-0.15)\times 10^{-3}\leavevmode\nobreak\ \mathrm{GeV}^{-2}$ for all LECs. With the current precision, our value for $C_{7}^{\mathrm{Wr}}$ in Eq. (57) is also compatible with the larger LMD+P estimate.

As mentioned above, the precision obtained for the normalization of the TFF from the lattice cannot compete with the PrimEx or even the PrimEx-II result. However, with our new estimate for $C_{7}^{\mathrm{Wr}}$ in Eq. (57), we can follow the same strategy as in Refs. Moussallam:1994xp ; Ananthanarayan:2002kj ; Kampf:2009tk and use it together with the experimental decay width $\Gamma(\eta\to\gamma\gamma)=0.515(18)\leavevmode\nobreak\ \mathrm{keV}$ pdg:2018 to determine the other LEC $C_{8}^{\mathrm{Wr}}$ that appears in the expression for the decay $\pi^{0}\to\gamma\gamma$ at NLO in $\chi$PT with the result

[TABLE]

This result is almost identical to the estimate $C_{8}^{\mathrm{Wr}}=0.58(20)\times 10^{-3}\leavevmode\nobreak\ \mathrm{GeV}^{-2}$ from Ref. Kampf:2009tk that was obtained by setting $C_{7}^{\mathrm{Wr}}=0$. We took over their estimate of 30% uncertainty from neglected terms of order $m_{s}^{2}$ in the chiral expansion. Using the expression for the decay amplitude $\pi^{0}\to\gamma\gamma$ given in Ref. Kampf:2009tk , including strong isospin breaking effects from $m_{u}\neq m_{d}$ and electromagnetic corrections, these changes in the LECs lead to the prediction

[TABLE]

confirming the value $\Gamma(\pi^{0}\to\gamma\gamma)=8.09(11)\leavevmode\nobreak\ \mathrm{eV}$ obtained in Ref. Kampf:2009tk . Therefore the tension with the preliminary result of the PrimEx-II experiment Primex2Mainz remains.

Another phenomenologically interesting quantity is the slope of the form factor at the origin. It is defined through

[TABLE]

We obtain

[TABLE]

where the first error is statistical and the second error is the systematic error associated with the $z$-expansion. Both determinations are compatible with each other and with the result obtained in Ref. Hoferichter:2018dmo using a dispersive framework (within our convention, their result reads $b_{\pi^{0}}=1.73(5)\leavevmode\nobreak\ \mathrm{GeV}^{2}$). The model independent determination is affected by a relatively large error which could be reduced in the future by adding large-volume ensembles which probe the low-$Q^{2}$ region.

We also try the method of Canterbury approximants proposed in Ref. Masjuan:2017tvw in the context of the pion TFF. It is a generalization of the Padé approximants for bivariate functions and we refer the reader to Ref. Masjuan:2015lca for an introduction to the subject. A general Canterbury approximant has the form

[TABLE]

with the convention $b_{00}=1$ and the relations $a_{nm}=a_{mn}$, $b_{nm}=b_{mn}$ reflecting the Bose symmetry of the photons. We consider both the diagonal and the subdiagonal sequences of approximants $C_{N}^{N}$ and $C_{N+1}^{N}$. In particular, the sequence $C_{N+1}^{N}$ automatically satisfies the short-distance behavior $\sim 1/Q^{2}$ if one sets $b_{NN}=0$. The lowest order elements for both sequences are

[TABLE]

The main drawbacks of this method are, first, the absence of a proof of convergence for the pion TFF and, second, the rapid growth of the number of fit parameters. As a consequence, for large values of $M$, the denominator can become singular at large virtualities (the presence of spurious poles in the spacelike domain) where there are no lattice data. In the following, we perform a global fit of the TFF, over the whole kinematic range, assuming a linear dependence of the coefficients $a_{ij}$ and $b_{ij}$ on the dimensionless parameter $\widetilde{y}$ and a quadratic dependence in the lattice spacing $a$. Both sequences already give a good $\chi^{2}$ for $N=1$. Using higher-order approximants only increases the statistical error or leads to unstable fits. The results for the two sequences are summarized in Table 4. We use as our final estimate the average between the approximants $C_{2}^{1}$ and $C_{1}^{1}$ and we take half the difference between the two sequences as an estimate for the systematic error : $\alpha=0.261(12)_{\mathrm{stat}}(8)_{\mathrm{syst}}\leavevmode\nobreak\ \mathrm{GeV}^{-1}$. This value is compatible with our determination via the $z$-expansion with a similar statistical error, but a systematic error which is difficult to estimate.

IV.3 Asymptotic behavior of the pion transition form factor

When at least one of the photon virtualities is large, and assuming (collinear) factorization, the TFF can be written as a convolution integral BL_3_papers

[TABLE]

where $T_{H}(x,Q_{1}^{2},Q_{2}^{2},\mu^{2})$ is a hard scattering kernel, calculable in perturbative QCD, and $\varphi_{\pi}(x;\mu^{2})$ is the nonperturbative pion distribution amplitude (DA) of twist two (normalized to one). The latter is scale dependent but otherwise universal and therefore plays a major role in the study of hard exclusive processes. The factorization and renormalization scales are chosen to be equal $\mu_{R}=\mu_{F}=\mu$ and of order $\overline{Q}^{2}=Q_{1}^{2}+Q_{2}^{2}$. Here, $x$ and $1-x$ are, respectively, the quark and antiquark longitudinal momentum fractions of the meson’s total momentum. Neglecting isospin breaking effects, the pion distribution amplitude is symmetric under the interchange of $x$ and $\overline{x}=1-x$ such that $\varphi_{\pi}(x)=\varphi_{\pi}(1-x)$. At leading twist and leading order in perturbative QCD, the hard scattering kernel is given by BL_3_papers

[TABLE]

where we have introduced the asymmetry parameter $\omega$ defined through

[TABLE]

In this section, the parameter $\omega$ should not be confused with the parameter $\omega_{1}$ introduced in Eq. (3). At next-to-leading order in perturbative QCD, the hard scattering kernel has been computed for all values of $\omega$ in the $\overline{{\rm MS}}$ scheme in Refs. delAguila:1981nk ; Braaten:1982yp ; Kadantseva:1985kb and the result is

[TABLE]

where an explicit expression for $t(\overline{x},\overline{w})$ with $\overline{w}=(1-\omega)/2$ is given by Eq. (5.2) of Ref. Braaten:1982yp . Then, using $C_{F}=4/3$ and the symmetry properties of the pion DA, Eq. (65) reduces to

[TABLE]

The asymptotic expression of the twist two pion DA, $\varphi_{\pi}^{\rm as}(x)=6x(1-x)$, leads to the OPE and BL asymptotic predictions for the double and single off-shell form factors respectively BL_3_papers ; Efremov:1979qk . However, this asymptotic result is expected to hold only at very large virtualities. Higher twist corrections, discussed in Refs. Novikov:1983jt ; Khodjamirian:1997tk , are given by

[TABLE]

where the parameter $\delta^{2}=0.20(2)\leavevmode\nobreak\ \mathrm{GeV}^{2}$ has been evaluated in Ref. Novikov:1983jt using QCD sum rules and more recently in Ref. Bali:2018spj using lattice QCD.

IV.3.1 Double-virtual form factor : Symmetric case

When both photons share the same virtualities ($\omega=0$) one can see from Eq. (69) that the asymptotic value of the TFF is not sensitive to the shape of the pion DA (one can show that $t(\overline{x},1/2)=-3/2$ is also independent of the quark longitudinal momentum fraction $x$) and the result reads

[TABLE]

Checking the validity of Eq. (71) at asymptotically large $Q^{2}$ provides a strong test of perturbative QCD. Here we will ask whether the given functional form describes our lattice data at moderately large $Q^{2}$.

IV.3.2 Double-virtual form factor : General case

When both photons share different virtualities, the situation is more difficult because the result depends on the precise shape of the pion DA, which is largely unknown, although a recent lattice calculation of the two lowest Gegenbauer moments has been carried out on a similar set of lattice ensembles in Ref. Bali:2019dqc . However, as pointed out in Ref. Diehl:2001dg , this dependence is small for small values of the asymmetry parameter $\omega$. At leading twist, we have

[TABLE]

where $a_{2}$ is the second coefficient in the expansion of the pion distribution amplitude in terms of Gegenbauer polynomials. For $Q_{2}^{2}=2Q_{1}^{2}$, where $\omega=1/3$, corrections to the asymptotic DA are of the order of $2\%$. We can therefore fit our lattice data in a much wider range without being sensitive to the actual shape of the pion DA.

IV.3.3 Fits

We perform a global fit using Eq. (71) where $\delta^{2}$ is considered as a free fit parameter and is allowed to vary linearly with $\widetilde{y}$ and $a^{2}$. We restrict the fit to virtualities with $\omega<1/3$ and $Q^{2}>Q^{2}_{\rm min}$ and use the four-loop running strong coupling in the $\overline{\rm MS}$ scheme vanRitbergen:1997va . Of course, one might question the applicability of perturbative QCD down to such low values of $|Q_{\rm min}|=1.2\leavevmode\nobreak\ \mathrm{GeV}$, even in Euclidean space, see Ref. Eidelman:1998vc for an analysis of nonperturbative effects in the Adler function. The results at the physical point and for different values of $Q^{2}_{\rm min}$ are listed in Table 5 and we quote as our final estimate the value at $Q^{2}_{\rm min}=2\leavevmode\nobreak\ \mathrm{GeV}^{2}$

[TABLE]

This result is compatible with the sum rule determination Novikov:1983jt and the lattice determination Bali:2018spj although with larger error.

IV.4 Determination of the pion-pole contribution to HLbL scattering in the $(g-2)_{\mu}$

The hadronic light-by-light scattering contribution to the anomalous magnetic moment of the muon is one of two dominant sources of uncertainty along with the hadronic vacuum polarization. For the HLbL scattering, a model independent dispersive approach has been proposed recently HLbL_DR and the dominant contribution is expected to be the pion-pole contribution which requires the pion TFF as input. Until recently, most estimates were based on model calculations where errors are difficult to estimate. In our previous work Gerardin:2016cqj , we provided the first lattice QCD calculation of the pion-pole contribution, with a statistical precision of about 12%. Recently, a data-driven determination based on the calculation of the TFF using dispersion theory was published Hoferichter:2018dmo . The result is compatible with our previous determination but with reduced uncertainty. In this section, we propose to use our results to improve our estimate of the pion-pole contribution from lattice QCD.

Following Jegerlehner:2009ry , the pion-pole contribution to the hadronic light-by-light scattering in the muon $(g-2)$ can be expressed as a three-dimensional integral involving two model-independent weight functions $w_{1}$ and $w_{2}$ and the product of a single-virtual TFF times a double-virtual TFF for arbitrary spacelike virtualities,

[TABLE]

The integrals run over the lengths $Q_{i}=|(Q_{i})_{\mu}|$, $i=1,2$ of the two Euclidean four-momentum vectors and the angle $\theta$ between them, $Q_{1}\cdot Q_{2}=Q_{1}Q_{2}\cos\theta$, where we defined $\tau=\cos\theta$. The analytical expressions for the dimensionless weight functions $w_{i}(Q_{1},Q_{2},\tau)$, $i=1,2$ can be found in Ref. Jegerlehner:2009ry . In particular, it was shown in Ref. Nyffeler:2016gnb that the relevant integration range involves spacelike virtualities below $2\leavevmode\nobreak\ \mathrm{GeV}^{2}$, precisely the kinematical region where we have lattice data.

IV.4.1 Phenomenological models

A first estimate of the pion-pole contribution is obtained using the LMD+V model, which provides a good description of our lattice data. The parameters at the physical point are determined from the global fit procedure described in Sec. IV.1 and we obtain

[TABLE]

where the error is statistical but includes the error from the continuum and chiral extrapolations. It can be compared with our previous estimate, $a_{\mu;\mathrm{LMD}+\mathrm{V}}^{\mathrm{HLbL};\pi^{0}}=(65.0\pm 8.3)\times 10^{-11}$, obtained with two dynamical quarks Gerardin:2016cqj . The slightly lower central value in Eq. (75) arises due to the slightly lower values of the parameters $\alpha$ and $\widetilde{h}_{2}$ emerging from the fit, but the results are in agreement within the quoted uncertainties. We also checked that computing $a_{\mu;\mathrm{LMD}+\mathrm{V}}^{\mathrm{HLbL};\pi^{0}}$ on each ensemble separately and then, in a second step, extrapolating the results to the continuum limit and at the physical pion mass using a fit linear in $\widetilde{y}$ and $a^{2}$ gives a similar result. In this case, we would obtain $(58.4\pm 2.4)\times 10^{-11}$. Since we are using a phenomenological model to describe the lattice data, the result could be biased and the error underestimated. Moreover, for the LMD+V model, we do not fit all the model parameters but made some assumptions on the second vector-resonance mass $M_{V_{2}}$ to stabilize the fit.

IV.4.2 Canterbury approximants

A second possibility is to use the results obtained from the method of the Canterbury approximants presented in Sec. IV.2. The short-distance constraints in Eqs (39) and (40) are explicitly implemented using the sequence with $M=N+1$ with the additional constraint $b_{NN}=0$. Using the approximant $(N,M)=(1,2)$, which provides a good description of our lattice data, we find at the physical point

[TABLE]

The next approximant leads to unstable fits. We note that this result is compatible with the LMD+V model determinations, which suggests that the model dependence of the central value is small. As for the LMD+V model, it is, however, difficult to assess a systematic error, especially because we are limited to rather low values of $M$.

IV.4.3 Final result : $z$-expansion

Finally, the $z$-expansion provides a systematically improvable determination. The choice for the function $P(Q_{1}^{2},Q_{2}^{2})$ in Eq. (41) guarantees the $1/Q^{2}$ falloff of the TFF in all directions, and, using the results of Sec. III.3 with $N=3$, we quote as our final result

[TABLE]

where the first error is statistical, the second is the systematic error inferred from the study in Appendix B and the third error comes from the disconnected contribution. In particular, the first error includes the error on the lattice spacing, the renormalization of the vector currents and the extrapolation to the physical point. In the last step of Eq. (77), we have added the errors in quadrature. Since the TFF has units of $\mathrm{GeV}^{-1}$, the relative scale-setting uncertainty on the pion-pole contribution is 2 times the relative error on the lattice spacing expressed in fm. The scale is known with a precision of $1\%$ which translates in an uncertainty of about 2% on $a_{\mu}^{\mathrm{HLbL};\pi^{0}}$. The uncertainty on the renormalization of the vector current is negligible. Therefore, the statistical precision of the correlation functions and the extrapolation to the physical point are the dominant sources of uncertainties in our calculation. They could be improved in a future calculation, in particular by including a lattice ensemble at the physical pion mass in the analysis. In contrast to phenomenological analyses, our result for $a_{\mu}^{\mathrm{HLbL};\pi^{0}}$ is significantly more accurate than our determination of ${\cal F}_{\pi^{0}\gamma^{*}\gamma^{*}}(0,0)$.

In the left panel of Fig. 13, we show the difference $\Delta a_{\mu}^{\mathrm{HLbL};\pi^{0};\mathrm{disc}}$ between the results obtained with and without including the disconnected contribution. At the SU(3)flavor symmetric point, the difference vanishes exactly, while it turns negative as the pion mass approaches its physical value. At a constant value of the trace of the quark-mass matrix, the disconnected contribution is proportional to $m_{s}-m_{l}$ close to the SU(3)flavor symmetric point. In practice, we parametrize it as being linear in $m_{K}^{2}-m_{\pi}^{2}$ and obtain $\Delta a_{\mu}^{\mathrm{HLbL};\pi^{0};\mathrm{disc}}=-1.0(0.3)\times 10^{-11}$ at the physical point. To be conservative, and because our data do not allow us to perform a precise continuum extrapolation, we associate 50% uncertainty to this contribution.

Similar to the case of the LMD+V model, we could perform the $z$-expansion on each ensemble separately and then extrapolate $a_{\mu}^{\mathrm{HLbL};\pi^{0}}$ to the physical point. We would obtain $(58.5\pm 4.0)\times 10^{-11}$, compatible with the results from the global fit. We note that the continuum and chiral extrapolations are very mild, as can be seen in the right panel of Fig. 13. We prefer the global fit method which reduces the number of fit parameters and relies on the chiral extrapolation of the pion TFF itself instead of the pion-pole contribution.

Our final value given by Eq. (77) is compatible with our previous determination Gerardin:2016cqj using two dynamical quarks but with an improved accuracy. It is also in good agreement with the data-driven determination recently published in Ref. Hoferichter:2018dmo and based on dispersion theory, with a similar precision to results published in Ref. Masjuan:2017tvw , based on a fit to experimental data using Canterbury approximants.

IV.4.4 Combination of lattice and experimental data

The experimental precision on the normalization of the TFF, dominated by the PrimEx experiment Larin:2010kq , is still better than our lattice determination. It is thus possible to reduce the error on $a_{\mu}^{\mathrm{HLbL};\pi^{0}}$ by exploiting this experimental measurement. We thus performed a constrained linear fit: to propagate the error on the normalization of the TFF, the constraint is imposed on each jackknife sample with a Gaussian distribution reproducing the PrimEx result. We obtain

[TABLE]

Two comments on these results. First, this value is slightly higher than in Eq. (77). This can be explained by our lower value for the normalization of the TFF. Second, the statistical error has been reduced by a factor of 1.7 which can be attributed to the small experimental error on the decay width, emphasizing the importance of precise lattice data at low virtualities. Thus, to further improve the lattice determination, it would be interesting to add large volume ensembles, which more tightly constrain the normalization of the TFF. The CLS ensemble E250 Mohler:2017wnb , at the physical mass and with a volume of 6.2 fm would be valuable.

V Conclusion

We have computed the neutral pion transition form factor with $N_{f}=2+1$ dynamical quarks at the physical point. The results are summarized in Table 3 where we provide the coefficients of the modified $z$-expansion and the associated correlation matrix. Our results are in good agreement with experimental data from CELLO and CLEO in the single-virtual case. In the double-virtual case, where no experimental data exist yet, we can compare our result with the recent dispersive analysis Hoferichter:2018dmo . Finally, when both virtual photons carry large virtualities, our result reproduces the asymptotic prediction from perturbative QCD once one includes the first-order $\alpha_{s}$ correction and higher-twist corrections.

In Sec. IV, we first compared our result with popular phenomenological models, often used to estimate the pion-pole contribution to hadronic light-by-light scattering in the muon $g-2$. The LMD+V model, which satisfies both the OPE and the Brodsky-Lepage constraints at short distances, provides a good parametrization over the whole kinematic range covered by our lattice data. The fit parameters turn out to be close to their phenomenological values obtained in Ref. Nyffeler:2009tw . However, our lattice data show significant deviations from the VMD or LMD models. Second, we have extracted the normalization of the TFF using either the phenomenological models or the more systematic parametrization based on the $z$-expansion. We reproduce the experimental result with a precision of 3.5% : $\alpha=0.264(9)\leavevmode\nobreak\ \mathrm{GeV}^{-1}$. This is an important benchmark of our calculation. The precision is not yet competitive with the experimental estimate from the PrimEx Collaboration Larin:2010kq , but we were able to extract the LEC appearing in the odd-intrinsic-parity sector of $\chi$PT at order $p^{6}$, $C_{7}^{\mathrm{Wr}}=0.16(18)\times 10^{-3}\leavevmode\nobreak\ \mathrm{GeV}^{-2}$. This value lies in-between some recent estimates based on various resonance Lagrangians that do or do not fulfill short-distance constraints from the OPE. Ref. Kampf:2009tk revisited the analysis of pion decay $\pi^{0}\to\gamma\gamma$ at NLO in $\chi$PT and assumed that this LEC vanishes. It turns out that with our new result from the lattice, the other relevant LEC at NLO does not change much compared to the analysis in Ref. Kampf:2009tk . We get $C_{8}^{\mathrm{Wr}}=0.56(17)\times 10^{-3}\leavevmode\nobreak\ \mathrm{GeV}^{-2}$ and essentially reproduce with $\Gamma(\pi^{0}\to\gamma\gamma)=8.07(10)\leavevmode\nobreak\ \mathrm{eV}$ the result given in that reference, using similar input from experiment and theory.

Finally, our model-independent lattice estimate for the pion-pole contribution to hadronic light-by-light scattering in $(g-2)_{\mu}$, given by Eq. (77), reads $a_{\mu}^{\mathrm{HLbL};\pi^{0}}=(59.7\pm 3.6)\times 10^{-11}$ and corresponds to a precision of 6%. This is our main result. The precision can be further improved by imposing the experimental constraint on the normalization of the TFF from the PrimEx experiment, and we obtain in this way the lattice and data-driven estimate $a_{\mu}^{\mathrm{HLbL};\pi^{0}}=(62.3\pm 2.3)\times 10^{-11}$ with a precision of 4%. This precision has already reached a sufficient level in view of the forthcoming experiment at Fermilab, for which a precision of $10\%$ on the theory estimate of the full HLbL contribution is desired.

In the future, we plan to use our result to estimate the dominant finite-size effects in the full lattice calculation of the HLbL contribution to the muon $g-2$. Since the QCD four-point correlation function will be computed on the same set of lattice ensembles, we should be able to constrain the tail of the integrand at large distances and thus reduce the statistical error. We also plan to include another lattice ensemble, with physical pion mass and large volume, to constrain even better the chiral and continuum extrapolations. The large volume should help us reach smaller virtualities and will constrain the normalization of the TFF even better. Finally, besides the dominant pion-pole contribution, it would be interesting to have a first principle estimate for the $\eta$ and $\eta^{\prime}$ pseudoscalar-pole contributions. According to model calculations Nyffeler:2016gnb , the size of these contributions is of the order of 20% of the pion-pole contribution and they are therefore not negligible. The same lattice methodology can be used to extract the relevant pseudoscalar transition form factors, even though disconnected diagrams will play a more important role. In that case, the weight functions appearing in Eq. (74) are peaked at slightly larger virtualities. However, since a precision of 20 or 30% should suffice, a lattice calculation should be feasible. Such a study would have a high impact since only sparse experimental data are available in the double-virtual case for spacelike momenta BaBar:2018zpn .

Acknowledgements.

We thank Martin Hoferichter and Bastian Kubis for helpful discussions. We are grateful for the access to the ensembles used here, made available to us through CLS, as well as to the samples of quark-disconnected loops generated by Konstantin Ottnad on a subset of these ensembles. The quark-connected correlation functions were computed on the platforms “Clover” at the Helmholtz-Institut Mainz and “Mogon II” at Johannes Gutenberg University Mainz. The CLS lattice ensembles and the quark-disconnected loops used here were partly generated through computing time provided by the Gauss Centre for Supercomputing (GCS) through the John von Neumann Institute for Computing (NIC) on the GCS share of the supercomputer JUQUEEN at Jülich Supercomputing Centre (JSC). This work was partly supported by the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme through Grant Agreement No. 771971-SIMDAMA, as well as by the Cluster of Excellence Precision Physics, Fundamental Interactions, and Structure of Matter (PRISMA+ EXC 2118/1) funded by the German Research Foundation (DFG) within the German Excellence Strategy (Project ID 39083149), as well as by the SFB 1044 funded by the DFG. A.N. is grateful to the Albert Einstein Center for Fundamental Physics (AEC) at the University of Bern for hospitality and support during the completion of this work.

Appendix A Analytic expressions of $\widetilde{A}_{\mu\nu}^{\rm VMD}(\tau)$ and $\widetilde{A}_{\mu\nu}^{\rm LMD}(\tau)$

In this appendix, we provide analytical expressions for the function $\widetilde{A}_{\mu\nu}^{\rm LMD}(\tau)$, introduced in Eq. (3), for a general pion momentum $\vec{p}$, assuming an LMD transition form factor Moussallam:1994xp ; Knecht:1999gb . The VMD case is simply obtained by setting $\beta=0$ in the equations below. Using the LMD model and Eq. (10), we obtain

[TABLE]

where $P^{E}_{\mu\nu}=i\epsilon_{\mu\nu 0i}p^{i}$ and $Q^{E}_{\mu\nu}=\epsilon_{\mu\nu i0}\,E_{\pi}\,q_{1}^{i}\,-i\epsilon_{\mu\nu ij}\,q_{1}^{i}\,p^{j}$ are independent of $\widetilde{\omega}$. The integrand has four distinct simple poles

[TABLE]

where we used the notations $k_{1}=\sqrt{M_{V}^{2}+|\vec{q}_{1}|^{2}}$ and $k_{2}=\sqrt{M_{V}^{2}+|\vec{q}_{2}|^{2}}$. Therefore

[TABLE]

Case $\tau>0$

[TABLE]

Case $\tau<0$

[TABLE]

The expression for $\widetilde{A}^{(1),\mathrm{LMD}}(\tau)$ is then obtained using Eq. (10). The function $\widetilde{A}_{\mu\nu}^{\mathrm{LMD}}(\tau)$ is continuous for all values of $\tau$ when $P^{E}_{\mu\nu}=0$ but is discontinuous at $\tau=0$ when $P^{E}_{\mu\nu}\neq 0$. However, the function has finite limits when $\tau\to 0^{\pm}$ respectively. A similar method can easily be applied to the LMD+V model.

Finally, we remark that Eqs. (82) and (83) indicate the type of intermediate states that contribute to the Euclidean three-point function $C_{\mu\nu}^{(3)}$ defined in Eq. (5), according to the LMD model. Taking into account $\vec{p}=\vec{q}_{1}+\vec{q}_{2}$ and the relation (4) between $\widetilde{A}_{\mu\nu}(\tau)$ and $C_{\mu\nu}^{(3)}$, we find that for $\tau>0$, either a vector meson with momentum $\vec{q}_{1}$ or a two-particle state consisting of a pion with momentum $\vec{p}$ and a vector meson with momentum $-\vec{q}_{2}$ is propagating between the two vector currents; for $\tau<0$, either a vector meson with momentum $\vec{q}_{2}$ or a two-particle state consisting of a pion with momentum $\vec{p}$ and a vector meson with momentum $-\vec{q}_{1}$ is propagating.

Appendix B $z$-expansion : systematic error

The transition form factor admits a modified double $z$-expansion given by Eq. (43) where the conformal variables $z_{1}$ and $z_{2}$ are given by

[TABLE]

As explained in the main text, $t_{0}$ is the optimal choice which reduces the maximum value of $|z_{k}|$ in the range $[0,Q_{\max}^{2}]$. Using $m_{\pi}=134.9\leavevmode\nobreak\ \mathrm{MeV}$ one finds

[TABLE]

well below one. The expansion parameter $z_{k}$ vanishes for $Q_{k}^{2}=-t_{0}$ (with $Q_{\max}^{2}=4\leavevmode\nobreak\ \mathrm{GeV}^{2}$, it corresponds to $t_{0}=-0.47\leavevmode\nobreak\ \mathrm{GeV}^{2}$) and is maximum at $z_{k}=0$ and $z_{k}=Q^{2}_{\max}$. Finally, the virtualities can be expressed in terms of the variables $z_{k}$

[TABLE]

In Sec. III.3, the coefficients of the double $z$-expansion were obtained by fitting the lattice data for a given value of $N$. The latter should be large enough to minimize the systematic error coming from the truncation of the sum but its value is also limited by the statistical precision of the data. We would like to estimate the systematic error induced by this truncation. Therefore, as a test, we fit the LMD+V TFF, defined using Eq. (48c) and the fit parameters in Eq. (53), with the $z$-expansion given by Eq. (43). The model satisfies both short-distance constraints given by Eqs. (39) and (40). For this purpose, we use a regular grid with $0\leq Q_{k}^{2}\leq Q_{\max}^{2}$ and a step size $\delta Q_{k}^{2}=0.025\leavevmode\nobreak\ \mathrm{GeV}^{2}$. In Table 6, and for different values of $N$ and $Q_{\max}^{2}$, we provide the results of the fit for the normalization of the TFF $\alpha$, the slope at the origin $b_{\pi}$ and the pion-pole contribution $a_{\mu}^{\mathrm{HLbL};\pi^{0}}$ which can be compared to the exact known results. We also provide the maximum deviation between the exact TFF and the fit. The level of agreement is illustrated in Fig. 14.

We conclude that using $N=3$ and $Q^{2}_{\max}=4\leavevmode\nobreak\ \mathrm{GeV}^{2}$ is already sufficient to get a precision below 1% for the TFF in the range $[0,Q^{2}_{\max}]$. The normalization of the TFF and the pion-pole contribution, $a_{\mu}^{\mathrm{HLbL};\pi^{0}}$, are also recovered within a precision below 1%. The slope parameter is obtained with a precision of 6%. We point out that, using the optimal value $t_{0}$, the parameters $z_{k}$ reach their maximal value at $Q_{k}^{2}=0$, precisely where the normalization of the TFF is obtained. Therefore, fitting the data in a wide range of virtualities is not the best method to determine the normalization of the TFF. However, the minimal value of $z_{k}$ is obtained at $Q^{2}_{k}=-t_{0}\approx 0.5\leavevmode\nobreak\ \mathrm{GeV}^{2}$, an optimal choice for the pion-pole contribution since the integrand is peaked around this value.

Appendix C Results of the $z$-expansion for some individual lattice ensembles

In Table 7, we provide the coefficients of the $z$-expansion for three ensembles used in this work: N101, N200 and D200. The correlation matrices are respectively given by Eqs (87), (88) and (89). The ensembles N200 and D200 have the same lattice spacing but different pion masses whereas N200 and N101 have similar pion masses but different lattice spacings. We stress that our main result in Sec. III.3 is based on a global fit, where all the ensembles are fitted simultaneously using the ansatz (44).

[TABLE]

[TABLE]

[TABLE]