Pseudoscalar-pole contribution to the $(g_{\mu}-2)$: a rational approach

TL;DR

This paper introduces a systematic, data-driven rational approximant method to accurately calculate the pseudoscalar-pole contribution to the muon's anomalous magnetic moment, reducing theoretical uncertainties for upcoming experiments.

Contribution

It presents a novel, model-independent rational approximation framework for calculating the pseudoscalar-pole contribution to the muon g-2, with a systematic uncertainty estimate based on extensive data.

Findings

Calculated $a_{ extrm{mu}}^{ extrm{HLbL};P}$ as 94.3(5.3) x 10^{-11}

Method achieves uncertainty below future experimental precision

Provides a systematic, data-based approach with uncertainty quantification

Abstract

We employ a mathematical framework based on rational approximants in order to calculate the pseudoscalar-pole piece of the hadronic light-by-light contribution to the anomalous magnetic moment of the muon, $a_{\mu}^{\textrm{HLbL};P}$. The method is systematic and data based, profiting from over 13 different collaborations, and able to ascribe, for the first time, a systematic uncertainty which provides for the model independence. As a result, we obtain $a_{\mu}^{\textrm{HLbL};P}=94.3(5.3) \times10^{-11}$, which uncertainty is well below the one foreseen at future experiments measuring the $(g_{\mu}-2)$.

Figures

Tables

Table 1. Table 1: The results for a μ HLbL ; π 0 × 10 11 superscript subscript 𝑎 𝜇 HLbL superscript 𝜋 0 superscript 10 11 a_{\mu}^{\textrm{HLbL};\pi^{0}}\times 10^{11} using the Regge and logarithmic models (last row) are compared to their C N + 1 N ​ ( Q 1 2 , Q 2 2 ) subscript superscript 𝐶 𝑁 𝑁 1 superscript subscript 𝑄 1 2 superscript subscript 𝑄 2 2 C^{N}_{N+1}(Q_{1}^{2},Q_{2}^{2}) sequence of approximants’ results. The LE row uses a pure low-energy reconstruction, whereas the OPE n rows incorporate high-energy constraints. The Fact row serves as an illustration of what a factorization approach would have yield. Finally, Fit OPE OPE {}^{\textrm{OPE}} row shows what a C N + 1 N ​ ( Q 1 2 , Q 2 2 ) subscript superscript 𝐶 𝑁 𝑁 1 superscript subscript 𝑄 1 2 superscript subscript 𝑄 2 2 C^{N}_{N+1}(Q_{1}^{2},Q_{2}^{2}) -like fitting function with the appropriate OPE behavior would lead. More details in the main text.

	Regge Model				Log Model
	$C_1^0$	$C_2^1$	$C_3^2$	$C_4^3$	$C_1^0$	$C_2^1$	$C_3^2$	$C_4^3$
LE	$55.2$	$59.7$	$60.4$	$60.6$	$56.7$	$64.4$	$66.1$	$66.8$
OPE0	$65.7$	$60.8$	$60.7$	$60.7$	$65.7$	$67.3$	$67.5$	$67.6$
OPE1	$-$	$60.6$	$60.7$	$60.7$	$65.7$	$67.3$	$67.5$	$67.6$
OPE2	$-$	$60.8$	$60.7$	$60.7$	$65.7$	$67.3$	$67.5$	$67.6$
Fact	$54.6$	$57.3$	$57.4$	$57.5$	$54.6$	$60.3$	$61.3$	$61.6$
Fit${}^{\text{OPE}}$	$66.3$	$62.7$	$61.1$	$60.8$	$79.6$	$71.9$	$69.3$	$68.4$
Exact	$60.7$				$67.6$

Table 2. Table 2: Our a μ HLbL ; P superscript subscript 𝑎 𝜇 HLbL 𝑃 a_{\mu}^{\textrm{HLbL};P} results in units of 10 − 11 superscript 10 11 10^{-11} . The second column is our simplest C 1 0 ​ ( Q 1 2 , Q 2 2 ) subscript superscript 𝐶 0 1 superscript subscript 𝑄 1 2 superscript subscript 𝑄 2 2 C^{0}_{1}(Q_{1}^{2},Q_{2}^{2}) approximant; the third and fourth columns refer to the C 2 1 ​ ( Q 1 2 , Q 2 2 ) subscript superscript 𝐶 1 2 superscript subscript 𝑄 1 2 superscript subscript 𝑄 2 2 C^{1}_{2}(Q_{1}^{2},Q_{2}^{2}) one, and stand for the lower and upper a P ; 1 , 1 subscript 𝑎 𝑃 1 1 a_{P;1,1} values respectively. See description in the text.

$a_{\mu }^{\text{HLbL;}P}$	$C_1^0$	$C_2^1\left[a_{P;1,1}^{\text{min}}\right]$	$C_2^1\left[a_{P;1,1}^{\text{max}}\right]$
${\pi }^0$	$65.3\left(1.4\right)\left(2.4\right)\left[2.8\right]$	$64.1\left(1.3\right)\left[1.3\right]$	$63.0\left(1.1\right)\left(0.5\right)\left[1.2\right]$
$\eta$	$17.1\left(0.6\right)\left(0.2\right)\left[0.6\right]$	$16.3\left(0.8\right)\left[0.8\right]$	$16.2\left(0.8\right)\left(0.6\right)\left[1.0\right]$
${\eta }^{\prime }$	$16.0\left(0.5\right)\left(0.3\right)\left[0.6\right]$	$14.7\left(0.7\right)\left[0.7\right]$	$14.3\left(0.5\right)\left(0.5\right)\left[0.7\right]$
Total	$98.4\left[2.9\right]$	$95.1\left[1.7\right]$	$93.5\left[1.7\right]$

Table 3. Table 3: The a μ HLbL ; π 0 superscript subscript 𝑎 𝜇 HLbL superscript 𝜋 0 a_{\mu}^{\textrm{HLbL};\pi^{0}} result ( 10 − 11 superscript 10 11 10^{-11} units) from different resonance-like approaches employed to approximate the Regge model which include up to n 𝑛 n resonances in Eq. ( 30 ). The exact result to compare with is represented in the last row. Find details in the text.

$n$	$1$	$2$	$3$	$4$
Res	$38.1$	$47.1$	$50.8$	$52.8$
Norm	$50.8$	$57.0$	$58.4$	$59.1$
Der	$50.8$	$57.8$	$59.4$	$59.9$
Fit	$55.9$	$67.2$	$58.3$	$65.4$
Fit${}^{\text{OPE}}$	$-$	$63.3$	$58.0$	$61.5$
Exact	$60.7$

Table 4. Table 4: The analog results to those for a μ HLbL; ​ π 0 superscript subscript 𝑎 𝜇 HLbL; superscript 𝜋 0 a_{\mu}^{\textrm{HLbL;}\pi^{0}} ( 10 − 11 superscript 10 11 10^{-11} units) in Table  2 , but employing the light-quark TFF. See details in the text.

$a_{\mu }^{\text{HLbL;}l.q.}$	$C_1^0$	$C_2^1\left[a_{P;1,1}^{\text{min}}\right]$	$C_2^1\left[a_{P;1,1}^{\text{max}}\right]$
$l.q.$	$60.4\left(1.5\right)\left(0.5\right)\left[1.6\right]$	$57.2\left(1.8\right)\left[1.8\right]$	$57.3\left(1.4\right)\left(1.0\right)\left[1.8\right]$
$l.q.norm$	$67.4\left(1.7\right)\left(0.5\right)\left[1.8\right]$	$63.8\left(2.0\right)\left[2.0\right]$	$63.9\left(1.6\right)\left(1.1\right)\left[1.9\right]$

Table 5. Table 5: The results for a μ HLbL; ​ P superscript subscript 𝑎 𝜇 HLbL; 𝑃 a_{\mu}^{\textrm{HLbL;}P} in units of 10 − 11 superscript 10 11 10^{-11} according to the procedure in Ref.  Melnikov and Vainshtein ( 2004 ) . The errors and labeling are identical to those in Table  2 .

$a_{\mu }^{\text{HLbL;}P}$	$C_1^0$	$C_2^1\left[a_{P;1,1}^{\text{min}}\right]$	$C_2^1\left[a_{P;1,1}^{\text{max}}\right]$
${\pi }^0$	$84.9\left(1.8\right)\left(2.6\right)\left[3.2\right]$	$82.8\left(1.7\right)\left[1.7\right]$	$80.9\left(1.3\right)\left(0.5\right)\left[1.4\right]$
$\eta$	$29.1\left(1.0\right)\left(0.3\right)\left[1.0\right]$	$27.3\left(1.4\right)\left[1.4\right]$	$26.9\left(1.5\right)\left(1.0\right)\left[1.8\right]$
${\eta }^{\prime }$	$30.4\left(1.0\right)\left(0.5\right)\left[1.1\right]$	$26.8\left(1.1\right)\left[1.1\right]$	$25.8\left(0.7\right)\left(0.9\right)\left[1.1\right]$
Total	$144.4\left[3.5\right]$	$136.9\left[2.5\right]$	$133.6\left[2.5\right]$

Table 6. Table 6: The values for the different parameters employed in calculating a μ HLbL; ​ P superscript subscript 𝑎 𝜇 HLbL; 𝑃 a_{\mu}^{\textrm{HLbL;}P} (second column) and the reference where they were extracted from (third column). F P ​ γ ​ γ subscript 𝐹 𝑃 𝛾 𝛾 F_{P\gamma\gamma} , P ∞ subscript 𝑃 P_{\infty} and δ π subscript 𝛿 𝜋 \delta_{\pi} are expressed in GeV -1 , GeV and GeV -2 units, respectively (additional parameters are dimensionless). The employed and expected future data are given in the fourth and fifth column, respectively.

Input	Value	Refs.	Data	Future
$F_{\pi \gamma \gamma }$	$0.2724\left(29\right)$		Olive et al. (2014)	Gasparian (2016); Babusci et al. (2012)
$b_{\pi }$	$0.0321\left(19\right)$	Masjuan (2012)	Behrend et al. (1991); Gronberg et al. (1998); Aubert et al. (2009); Uehara et al. (2012); Adlarson et al. (2016b); Lazzeroni et al. (2016)	Babusci et al. (2012); Adl (2014); Gan (2015)
$c_{\pi }$	$0.00104\left(22\right)$	Masjuan (2012)	Behrend et al. (1991); Gronberg et al. (1998); Aubert et al. (2009); Uehara et al. (2012); Adlarson et al. (2016b); Lazzeroni et al. (2016)	Babusci et al. (2012); Adl (2014); Gan (2015)
$d_{\pi }$				Babusci et al. (2012); Adl (2014); Gan (2015)
$a_{\pi ;1,1}$				Adl (2014)
${\pi }_{\infty }$	$2F_{\pi }$	Olive et al. (2014)
${\delta }_{\pi }^2$	$0.20\left(2\right)$	Novikov et al. (1984); Jegerlehner and Nyffeler (2009)
$F_{\eta \gamma \gamma }$	$0.2738\left(47\right)$		Olive et al. (2014)	Adl (2014)
$b_{\eta }$	$0.572\left(8\right)$	Escribano and Gonzalez-Solis (2015); Escribano et al. (2014)	Behrend et al. (1991); Gronberg et al. (1998); del Amo Sanchez et al. (2011); Arnaldi et al. (2009); Berghauser et al. (2011); Aguar-Bartolome et al. (2014); Arnaldi et al. (2016); Adlarson et al. (2016a)	Adl (2014); Gan (2015, 2015)
$c_{\eta }$	$0.333\left(9\right)$	Escribano and Gonzalez-Solis (2015); Escribano et al. (2014)	Behrend et al. (1991); Gronberg et al. (1998); del Amo Sanchez et al. (2011); Arnaldi et al. (2009); Berghauser et al. (2011); Aguar-Bartolome et al. (2014); Arnaldi et al. (2016); Adlarson et al. (2016a)	Adl (2014); Gan (2015)
$d_{\eta }$	$0.195\left(20\right)$	Escribano and Gonzalez-Solis (2015); Escribano et al. (2014)	Behrend et al. (1991); Gronberg et al. (1998); del Amo Sanchez et al. (2011); Arnaldi et al. (2009); Berghauser et al. (2011); Aguar-Bartolome et al. (2014); Arnaldi et al. (2016); Adlarson et al. (2016a)	Adl (2014); Gan (2015)
${\eta }_{\infty }$	$0.180\left(12\right)$	Escribano and Gonzalez-Solis (2015); Escribano et al. (2014)	Behrend et al. (1991); Gronberg et al. (1998); del Amo Sanchez et al. (2011); Arnaldi et al. (2009); Berghauser et al. (2011); Aguar-Bartolome et al. (2014); Arnaldi et al. (2016); Adlarson et al. (2016a)	Adl (2014); Gan (2015)
$F_{{\eta }^{\prime }\gamma \gamma }$	$0.3437\left(55\right)$		Olive et al. (2014)	Dudek et al. (2012); Adl (2014)
$b_{{\eta }^{\prime }}$	$1.31\left(3\right)$	Escribano et al. (2016, 2014)	Behrend et al. (1991); Gronberg et al. (1998); Acciarri et al. (1998); del Amo Sanchez et al. (2011); Ablikim et al. (2015)	Dudek et al. (2012); Adl (2014)
$c_{{\eta }^{\prime }}$	$1.74\left(9\right)$	Escribano et al. (2016, 2014)	Behrend et al. (1991); Gronberg et al. (1998); Acciarri et al. (1998); del Amo Sanchez et al. (2011); Ablikim et al. (2015)	Dudek et al. (2012); Adl (2014)
$d_{{\eta }^{\prime }}$	$2.30\left(22\right)$	Escribano et al. (2016, 2014)	Behrend et al. (1991); Gronberg et al. (1998); Acciarri et al. (1998); del Amo Sanchez et al. (2011); Ablikim et al. (2015)	Dudek et al. (2012); Adl (2014)
${\eta }_{\infty }^{\prime }$	$0.255\left(4\right)$	Escribano et al. (2016, 2014)	Behrend et al. (1991); Gronberg et al. (1998); Acciarri et al. (1998); del Amo Sanchez et al. (2011); Ablikim et al. (2015)	Dudek et al. (2012); Adl (2014)

Table 7. Table 7: Parameters of the TFF parameterization Eq. ( 35 ) obatined from our fitting procedure.

$t_0$	$t_1$	$t_2$	$r_1$	$r_2$	$r_3$
$0.276$	$-0.00024$	$0.00024$	$1.729$	$-0.0285$	$0.0013$

Equations

$${\Pi }_{\mu \nu \lambda \rho }(q_1,q_2,q_3)$$

$$\times ⟨\Omega |T\{j_{\mu }(x)j_{\nu }(y)j_{\lambda }(z)j_{\rho }(0)\}|\Omega ⟩,$$

$$a_{\mu }^{\text{HLbL}}=$$

$$\displaystyle(\not{p}+m_{\mu})\gamma^{\mu}(\not{p}-\not{q}_{1}\!+m_{\mu})\gamma^{\nu}(\not{p}-\not{q}_{1}\!-\not{q}_{2}\!+m_{\mu})\gamma^{\lambda}\Big{)}\!\left[\frac{\partial}{\partial k^{\rho}}\Pi_{\mu\nu\lambda\sigma}(q_{1},q_{2},k-q_{1}-q_{2})\right]_{k\rightarrow 0}.$$

$$\begin{gathered} {\Pi }^{\mu \nu \lambda \rho }(q_1,q_2,q_3)=\int d^{4}x\int d^{4}z{e}^{i{q}_1\cdot x}{e}^{i{q}_3\cdot z}\frac{i}{q^2-m_P^2+iϵ} \\ \times ⟨0|T\{j^{\mu }(x)j^{\nu }(0)\}|P⟩⟨P|T\{j^{\lambda }(z)j^{\rho }(0)\}|0⟩+\text{OT}, \end{gathered}$$

$$i{\mathcal{M}}_{P\to {\gamma }^{\ast }{\gamma }^{\ast }}^{\mu \nu }$$

$$\equiv -i{ϵ}^{\mu \nu \rho \sigma }{q}_{1\rho }{q}_{2\sigma }{F}_{P{\gamma }^{\ast }{\gamma }^{\ast }}(q_1^2,q_2^2),$$

$$a_{\mu }^{\text{HLbL};P}=$$

$$\times \left[\frac{F_1{I}_1(Q_1,Q_2,t)}{Q_2^2+m_P^2}+\frac{F_2{I}_2(Q_1,Q_2,t)}{Q_3^2+m_P^2}\right]$$

$$\equiv$$

$$F_1$$

$$F_2$$

$$w_2(Q_1,Q_2,t)\sim \frac{Q_{1(2)}^2}{m_P^2}\left(a(t)+b(t)\frac{Q_{2(1)}}{m_{\mu }}+...\right)+...,$$

$$\underset{Q_{1(2)}^2\to \infty }{\lim }{w}_2(Q_1,Q_2,t)$$

$$\underset{Q^2\to \infty }{\lim }{w}_2(Q,Q,t)$$

$$\underset{Q_1^2\to \infty }{\lim }{w}_1(Q_1,Q_2,t)\sim Q_1^{-1},$$

$$\underset{Q_2^2\to \infty }{\lim }{w}_1(Q_1,Q_2,t)\sim Q_2^{-2},$$

$$\underset{Q^2\to \infty }{\lim }{w}_1(Q,Q,t)\sim Q^{-2}$$

$$f(x,y)=\sum _{i,j}{c}_{i,j}{x}^i{y}^j(c_{i,j}=c_{j,i}),$$

$$C_M^N(x,y)=\frac{R_N(x,y)}{Q_M(x,y)}=\frac{{\sum }_{i,j=0}^N{a}_{i,j}{x}^i{y}^j}{{\sum }_{i,j=0}^M{b}_{i,j}{x}^i{y}^j}$$

$$\sum _{i=0}^{\alpha }\sum _{j=0}^{\beta }{b}_{i,j}{c}_{\alpha -i,\beta -j}=a_{\alpha ,\beta }\text{for }(\alpha ,\beta )\in \mathcal{N},$$

$$\sum _{i=0}^{\text{min}(\alpha ,M)}\sum _{j=0}^{\text{min}(\beta ,M)}{b}_{i,j}{c}_{\alpha -i,\beta -j}=a_{\alpha ,\beta }\text{for }(\alpha ,\beta )\in \mathcal{E}\notin \mathcal{N},$$

$$F_{P{\gamma }^{\ast }{\gamma }^{\ast }}^{\text{Regge(Log)}}(Q^2,Q^2)={\mathcal{C}}_1{Q}^{-2}+{\mathcal{C}}_2{Q}^{-4}+....$$

$$\begin{gathered} F_{P{\gamma }^{\ast }{\gamma }^{\ast }}(Q_1^2,Q_2^2)=F_{P\gamma \gamma }(0,0)(1-\frac{b_P}{m_P^2}(Q_1^2+Q_2^2) \\ +\frac{c_P}{m_P^4}(Q_1^4+Q_2^4)+\frac{a_{P;1,1}}{m_P^4}{Q}_1^2{Q}_2^2+...). \end{gathered}$$

$$F_{P\gamma \gamma }(0,0)=\frac{N_c}{4{\pi }^{2}F}\mathrm{tr}({\mathcal{Q}}^2{\lambda }_P).$$

$$\underset{Q^2\to \infty }{\lim }{F}_{P{\gamma }^{\ast }\gamma }(Q^2,0)=P_{\infty }{Q}^{-2}+\mathcal{O}(Q^{-4}),$$

$$\underset{Q^2\to \infty }{\lim }{F}_{P{\gamma }^{\ast }{\gamma }^{\ast }}(Q^2,Q^2)=\frac{P_{\infty }}{3}\left(\frac{1}{Q^2}-\frac{8}{9}\frac{{\delta }_P^2}{Q^4}\right)+\mathcal{O}(Q^{-6}),$$

$$C_1^0(Q_1^2,Q_2^2)$$

$$C_2^1(Q_1^2,Q_2^2)$$

$$a_{\mu }^{\text{HLbL};P}$$

$$\to 94.3(5.3)\times 10^{-11},$$

$$\begin{gathered} {\Pi }_{\mu \nu \lambda \rho }^{P-\text{pole}}(q_1,q_2,q_3)= \\ i\frac{F_{P{\gamma }^{\ast }{\gamma }^{\ast }}(q_1^2,q_2^2)F_{P{\gamma }^{\ast }{\gamma }^{\ast }}(q_3^2,k^2)}{(q_1+q_2{)}^2-m_P^2}{ϵ}_{\mu \nu \alpha \beta }{q}_1^{\alpha }{q}_2^{\beta }{ϵ}_{\lambda \rho \sigma \tau }{q}_3^{\sigma }{k}^{\tau } \\ +i\frac{F_{P{\gamma }^{\ast }{\gamma }^{\ast }}(q_1^2,k^2)F_{P{\gamma }^{\ast }{\gamma }^{\ast }}(q_3^2,q_2^2)}{(q_2+q_3{)}^2-m_P^2}{ϵ}_{\mu \rho \alpha \beta }{q}_1^{\alpha }{k}^{\beta }{ϵ}_{\nu \lambda \sigma \tau }{q}_2^{\sigma }{q}_3^{\tau } \\ +i\frac{F_{P{\gamma }^{\ast }{\gamma }^{\ast }}(q_1^2,q_3^2)F_{P{\gamma }^{\ast }{\gamma }^{\ast }}(k^2,q_2^2)}{(q_1+q_3{)}^2-m_P^2}{ϵ}_{\mu \lambda \alpha \beta }{q}_1^{\alpha }{q}_3^{\beta }{ϵ}_{\nu \rho \sigma \tau }{q}_2^{\sigma }{k}^{\tau }, \end{gathered}$$

$$I_{1}(Q_{1},Q_{2},t)=\frac{-1}{m_{\mu}^{2}Q_{3}^{2}}\Bigg{[}\frac{4m_{\mu}^{2}t}{Q_{1}Q_{2}}\\ +(1-R_{m_{1}})\bigg{(}\frac{2Q_{1}t}{Q_{2}}\-4(1-t^{2})\bigg{)}-(1-R_{m_{1}})^{2}\frac{Q_{1}t}{Q_{2}}\\ -8X(Q_{1},Q_{2},t)(Q_{2}^{2}-2m_{\mu}^{2})(1-t^{2})\Bigg{]},$$

$$I_{2}(Q_{1},Q_{2},t)=\frac{-1}{m_{\mu}^{2}Q_{3}^{2}}\Bigg{[}2(1-R_{m_{1}})\bigg{(}\frac{Q_{1}t}{Q_{2}}+1\bigg{)}\\ +2(1-R_{m_{2}})\bigg{(}\frac{Q_{2}t}{Q_{1}}+1\bigg{)}\\ +4X(Q_{1},Q_{2},t)\left(Q_{3}^{2}+2m_{\mu}^{2}(1-t^{2})\right)\Bigg{]},$$

$$X(Q_1,Q_2,t)=$$

$$z=$$

$$R_{m_i}=$$

Full text

I Introduction

The anomalous magnetic moment of the muon, $a_{\mu}\equiv(g_{\mu}-2)/2$, represents one of our finest tests of the Standard Model (SM) of particle physics and its most recent measurement, $a_{\mu}^{\textrm{exp}}=116592091(63)\times 10^{-11}$ Bennett et al. (2006); Olive et al. (2014), has reached the astonishing precision of $0.54$ ppm. At such precision $a_{\mu}$ does not only provide a beautiful test of our understanding of elementary interactions, but represents an interesting probe of physics beyond the SM. Indeed, there is at present a discrepancy among experiment and theory of around $3\sigma$ Jegerlehner and Nyffeler (2009)111The progress on this field is captured in at least three recent dedicated workshops on the $(g_{\mu}-2)$ Masjuan and Venanzoni (2013); Blum et al. (2014); D’Ambrosio et al. (2016). . For this reason, two new experiments have been projected both at Fermilab Lee Roberts (2011) and J-PARC Mibe (2010), which expect to measure $a_{\mu}$ at a precision of around $0.14$ ppm and would shed light on the nature of the present discrepancy. However, such a tremendous effort on reducing the experimental uncertainty would be in vain unless the current theoretical calculations would reach a similar accuracy. At present, the theoretical error is dominated by two different hadronic contributions: the hadronic vacuum polarization (HVP), which amounts to $58$ ppm to $a_{\mu}$ and with an uncertainty of around $0.36$ ppm Jegerlehner and Nyffeler (2009), and the hadronic light-by-light (HLbL), which amounts to $0.87$ ppm to $a_{\mu}$ and with an uncertainty of around $0.33$ ppm Jegerlehner and Nyffeler (2009), leading to a total theoretical error reading $0.49$ ppm. These calculations, involving complicated loop integrals, are hindered via the non-perturbative hadronic physics dominating the loop integrals. Fortunately, such complications can be overcome for the dominant contribution, the HVP, since it is related through the optical theorem to the $\sigma(e^{+}e^{-}\to\textrm{hadrons})$ cross section. It is expected that, in the near future, the ongoing experimental program will allow to reduce the HVP errors according to what future experiments require. Still, such effort would be fruitless unless a similar reduction in the precision of the HLbL is achieved, which would otherwise dominate the theoretical SM uncertainty and make the experimental efforts pointless.

By contrast to the HVP, the HLbL entails a much richer structure that avoids an easy connection to data. As a consequence, the existing calculations have required certain modeling and approximation procedures. Ascribing them a systematic error is a difficult task, but the variety of the present results Jegerlehner and Nyffeler (2009) already suggests an error which is potentially larger than future experiments’ precision, which demands a new, more accurate and less model-dependent evaluation. Among the different contributions to the hadronic-light-by-light, the pseudoscalar-pole seems to dominate the full quantity, requiring therefore the best precision. Fortunately, such quantity can be rigorously defined in a quantum field theory and related to the pseudoscalar transition form factors (TFFs), which are observable quantities. This offers an opportunity to perform a data-driven approach for this contribution provided that a reliable method is established. In this work, we discuss a novel method based on Canterbury approximants (bivariate Padé approximants) which provides a mathematically and data-based description for the involved TFFs in the space-like (SL) region, allowing then for a model-independent calculation for the pseudoscalar-pole contribution to the hadronic-light-by-light which includes, for the first time, the sought systematic error which has been missing in previous calculations. Special attention on the reasons which justify the use of Canterbury approximants in contrast to previous resonance approaches is given.

The paper is structured as follows: first, we briefly introduce the most general HLbL contribution to $a_{\mu}$ and its general features in Sect. II. Its main piece, the pseudoscalar pole, is presented in Sect. III; this includes also a brief overview of current theoretical approaches and motivates the reasons for our new study based on Canterbury approximants. These are subsequently presented in Sect. IV. Finally, we give our results in Sect. V and discuss the role of future data in Sect. VI. Much information is relegated to the appendices, including, among others, the impact of $P\to\bar{\ell}\ell$ decays, a discussion concerning the light-quark TFF, comments concerning the pseudoscalar-exchange approach, a discussion on dispersion relations, and our most-updated data input profitting from over 13 different collaborations.

II Hadronic light-by-light contribution

The HLbL contribution to $a_{\mu}$ (cf. Refs. Jegerlehner and Nyffeler (2009); Prades et al. (2009)) is depicted in Fig. 1, where the gray blob represents the HLbL tensor defined as

[TABLE]

with all the momenta, $q_{i}$, outgoing. The resulting contribution to $a_{\mu}$ can be expressed in terms of this hadronic quantity using projection techniques that allow to specialize to the kinematical limit relevant to $a_{\mu}$ $(k\to 0)$ in advance, obtaining Knecht and Nyffeler (2002):

[TABLE]

Such quantity requires as input the HLbL tensor at all energy scales, including the (non-perturbative) low energies, which turn out to play the major role in numerical calculations as we will illustrate for the particular case of the $\pi^{0}$-pole contribution. Qualitatively, this can be understood, after Wick rotation, from the propagators in Eq. (2).

A full description of the HLbL tensor is far from a trivial task. As an illustration, while the dominant HVP requires a scalar function depending on a single energy scale, the HLbL tensor involves 138 scalar functions222That number is however reduced after considering the Ward identities and the kinematic configuration relevant to $a_{\mu}$ Bijnens et al. (1996); Colangelo et al. (2015). and six scalar variables (the four photon virtualities, $\{q_{i}^{2}\}_{i=1-4}$, and 2 Mandelstam variables).

In order to deal with such object, E. de Rafael proposed more than twenty years ago de Rafael (1994) to split the most relevant contributions to this tensor according to a combined expansion in terms of the chiral and the large-number-of-colors ($N_{c}$) limits of QCD; whereas large-$N_{c}$ represents the only known perturbative approach to QCD, the chiral expansion, in powers of momentum ($q_{i}^{2}$), helps to identify those contributions which play the main role at low energies — the most important in the calculation. According to this proposal, the leading contributions to the HLbL tensor are the pseudoscalar loop contributions, at order $\mathcal{O}(N_{c}^{0},q^{4})$, and the pseudoscalar-pole terms, at order $\mathcal{O}(N_{c},q^{6})$. Following these ideas, subleading contributions will account for heavier resonances and the continuum quark-loop contributions, all of them of order $\mathcal{O}(N_{c},q^{8})$. It remains then the task to calculate all the relevant contributions as accurate and precise as possible.

Actually, most of the results in the literature follow de Rafael’s proposal (see Refs. Hayakawa and Kinoshita (1998); Bijnens et al. (2002); Knecht et al. (2002); Blokland et al. (2002); Melnikov and Vainshtein (2004); Dorokhov and Broniowski (2008); Nyffeler (2009); Jegerlehner and Nyffeler (2009); Prades et al. (2009); Hong and Kim (2009); Cappiello et al. (2011); Kampf and Novotny (2011); Masjuan (2012); Bijnens and Abyaneh (2012); Escribano et al. (2014); Roig et al. (2014); Dorokhov et al. (2015), including full and partial contributions to $a_{\mu}^{\textrm{HLbL}}$) finding values for $a_{\mu}^{\textrm{HLbL}}$ between basically $6\times 10^{-10}$ and up to almost $14\times 10^{-10}$. Among them, the Jegerlehner and Nyffeler review Jegerlehner and Nyffeler (2009), quoting $(11.4\pm 4.0)\times 10^{-10}$, and the Glasgow consensus Prades et al. (2009), written by Prades, de Rafael, and Vainshtein, and quoting $(10.5\pm 2.6)\times 10^{-10}$, represent in our opinion the standard reference values for the HLbL. They agree well since they only differ by few subtleties: they both used the model from Knecht and Nyffeler Knecht and Nyffeler (2002) to account for the dominant contribution, but differ on how to implement the high-energy QCD constrains and the error propagation. Neither of both approaches contain however systematic errors from the chiral and large-Nc limits Masjuan (2012); Masjuan and Vanderhaeghen (2015); Masjuan and Peris (2007, 2008), which are difficult to estimate (cf. also the discussion in Ref. Masjuan (2015)). An attempt to consider a dispersion relation to calculate a subset of the pieces appearing in de Rafael’s proposal has been recently proposed in Refs. Colangelo et al. (2014); Pauk and Vanderhaeghen (2014), but numerics are still to come. All in all, even though the QCD features for the HLbL are well understood Jegerlehner and Nyffeler (2009); Prades et al. (2009), the details of the particular calculations are important to get the numerical result to the final required precision.

Alternative approaches to calculate the HLbL contribution exist as well. As an example, a ballpark estimate can be obtained using the analytical result for the heavy quark loop contribution to the HLbL in the lines of Refs. Pivovarov (2003); Erler and Toledo Sanchez (2006); Boughezal and Melnikov (2011); Greynat and de Rafael (2012); Masjuan and Vanderhaeghen (2015), finding somewhat higher values around $(12\div 17)\times 10^{-10}$. Likewise, the Dyson-Schwinger approach adopted in Ref. Fischer et al. (2011); Goecke et al. (2011, 2013) and the lattice QCD simulations in Refs. Blum et al. (2015, 2016, 2017); Green et al. (2015) are also approaches which do not follow the same procedure. Whereas yet incomplete and with some progress still required, promising advances have been reported already Goecke et al. (2011); Eichmann et al. (2016); Blum et al. (2015, 2016, 2017).

Despite the combined chiral and large-$N_{c}$ counting yields the pseudoscalar loop contributions and the pseudoscalar-pole terms as the leading contributions to the HLbL tensor, suggesting a similar size, the situation is more subtle. The more careful discussion in Ref. Melnikov and Vainshtein (2004) observed that the typical size of the momentum running in the pion loop turns out to be of order $4m_{\pi}$ Melnikov and Vainshtein (2004), which implies a slow convergence of the chiral expansion with non-leading terms logarithmically enhanced. This reduces the pion loop contribution with respect to the pion pole. Phenomenologically, the former is found to be around four times smaller than the latter Jegerlehner and Nyffeler (2009), which becomes thereby the most relevant contribution to be calculated among the different HLbL contributions. Since such contribution is typically found to be of order $10^{-9}$, in order to meet the $0.14$ ppm precision of future experiments, a precision below $10\%$ is a priori desired, which is beyond traditional approaches’ performance333As we illustrate in the following sections, the current values Jegerlehner and Nyffeler (2009); Prades et al. (2009) entail potential large systematic errors above $10\%$..

In the following section we outline concisely what this contribution refers to and its relation to the pseudoscalar TFFs. In addition, the relevant kinematical regions of interest are identified, whereby the requirements that a TFF parameterization necessitates are obtained, which will naturally motivate our approach employing Canterbury approximants.

III Pseudoscalar-pole contribution

The pseudoscalar-pole contribution to the HLbL tensor can be easily obtained within the language of Green’s functions. Inserting the identity as a sum over the QCD spectrum ($\mathds{1}=\sum_{X}\int d\Pi_{X}\ket{X}\bra{X}$, with $\ket{X}$ on-shell intermediate hadronic states) within the HLbL tensor Eq. (1), it is obtained that such function exhibits well-isolated poles for the lightest pseudoscalar states $P=\{\pi^{0},\eta,\eta^{\prime}\}\in X$ which contribution to $a_{\mu}^{\textrm{HLbL}}$ is depicted in Fig. 2. This is, the HLbL tensor can be expressed as (sum over $P$ assumed)

[TABLE]

where $q=q_{1}+q_{2}=k-q_{3}$ and with OT referring both, to different time orderings (crossed $t$ and $u$ channels) and to additional contributions from the QCD spectrum ($X\neq P$ and not necessarily of pseudoscalar nature) which do not become singular as $q^{2}\to m_{P}^{2}$. Besides, for the on-shell pseudoscalar states, the matrix elements defined above are related to the pseudoscalar TFFs (gray blobs in Fig. 2), defined as444Note that the $(ie)^{2}$ coupling is already implicit in Eq. (2).

[TABLE]

with $\epsilon^{0123}=+1$. As a consequence, the contribution from the pseudoscalar poles can be calculated as model-independent as the employed TFF description is, without incurring in any ambiguity.

OT in Eq. (3) defines heavier states — which could be hardly included in this way given their widths Olive et al. (2014) — and continuum contributions. These processes cannot however induce a non-analytic behavior at energies close to the lightest pseudoscalar masses, which allow to disentangle the pseudoscalar-pole contribution, as previously said, unambiguously.

At this point, it is worth to make a brief digression and comment on what is referred as the pseudoscalar-exchange contribution, an analytically similar but theoretically different contribution to the pseudoscalar-pole one. The authors from Ref. Melnikov and Vainshtein (2004) realized that the pseudoscalar-pole contribution could not reproduce the high-energy QCD constraints imposed by the OPE for the full HLbL tensor. However, by setting the TFF involving the external photon to a constant one, such constraint can be satisfied. Their approach, as they mention Melnikov and Vainshtein (2004), was meant as a model to interpolate the full HLbL tensor from the low to the high energies, which includes, effectively, excited pseudoscalar resonances. Later on, Ref. Jegerlehner and Nyffeler (2009) criticized such approach and introduced what is known as the off-shell TFF. This approach intends to consider all the pseudoscalar contributions (i.e. heavier resonances and continuum) into an off-shell pseudoscalar TFF by connecting with the $\langle VVP\rangle$ Green’s function and imposing its well-known high-energy behavior Melnikov and Vainshtein (2004). We note that such procedure cannot be rigorously derived as a pole-contribution, and could be thought as a model interpolating the low- and high-energy behavior for the exchange of pseudoscalar-like resonances in Eq. (3).

It is not our intention to discuss how the high energies of the HLbL tensor should be implemented in terms of the pseudoscalar-pole contribution; our more modest concern is to discuss a model-independent data-based description for the pseudoscalar-pole contribution to the $a_{\mu}^{\textrm{HLbL}}$, inspired by new and forthcoming experimental results, lattice simulations, and dispersive representations.

The pseudoscalar-pole contribution to the HLbL tensor (Eq. (3)) is given explicitly in Eq. (23) and involves the $s,t$ and $u$ channels as illustrated in Fig. 2. Inserting this into Eq. (2), performing the Wick rotation and using the Gegenbauer technique to perform angular integrations Knecht and Nyffeler (2002), the pseudoscalar-pole contribution to $a_{\mu}^{\textrm{HLbL}}$ can be expressed as an integral over one angular variable and two space-like momenta Jegerlehner and Nyffeler (2009):

[TABLE]

where $F_{1}$ and $F_{2}$ are defined as

[TABLE]

with $Q_{3}^{2}=Q_{1}^{2}+Q_{2}^{2}+2Q_{1}Q_{2}t$. The $I_{1,2}(Q_{1},Q_{2},t)$ functions are defined in Appendix A and $w_{1,2}\equiv w_{1,2}(Q_{1},Q_{2},t)$. $F_{1,2}$ are products of the relevant TFFs appearing in the blobs in Fig. 2 and are to be evaluated in the SL region. Notice the product of a single- and a double-virtual TFF. The two integrands in Eq. (5) for the $\pi^{0}$ case are shown in Fig. 3 for a constant TFF ($F_{P\gamma^{*}\gamma^{*}}(Q_{1}^{2},Q_{2}^{2})\equiv 1$) and $t=0$ along with the $\eta^{\prime}$ case for the first integrand alone — similar results hold for different $t$ values555For a thorough description of these integrands, the interested reader is referred to Ref. Nyffeler (2016)..

In general, for a given pseudoscalar, the integrand involving $w_{2}$ is an order of magnitude smaller than that involving $w_{1}$ (cf. left vs. center panels in Fig. 3) and is peaked at low energies. Particularly, in such regime it can be described as

[TABLE]

which hints the large rising close to $Q_{i}^{2}\sim m_{P,\mu}^{2}$ (a similar behavior holds for $w_{1}$), whereas the high-energy behavior reads Nyffeler (2016)

[TABLE]

The integral involving $w_{2}$ remains therefore finite even for a constant TFF. As a consequence, only a precise description for the TFFs at very low SL energies is required. By contrast, the dominant integrand, which involves $w_{1}$, behaves as Nyffeler (2016)

[TABLE]

and constitutes therefore a divergent integral for a constant TFF. Moreover, as it can be observed from Fig. 3 (left and right panels), despite its peak at low energies, this integrand is sensitive to the region above $1$ GeV. This is specially important for heavier pseudoscalars such as the $\eta^{\prime}$ as it can be observed in Fig. 3 right. In these cases, the low-energy peak is less pronounced and the tail is relevant up to energies beyond $2$ GeV666The $\pi^{0}$ features a more pronounced peak as compared to the $\eta$ and $\eta^{\prime}$ and provides the main contribution. This is related to the chiral enhancement $\propto Q^{-2}$ from the pseudoscalar propagator, which is stronger for the $\pi^{0}$ given its mass.

Actually, for the TFFs we employ in Section V, the integral on $Q_{1}$ and $Q_{2}$ for $\pi^{0},\eta$ and $\eta^{\prime}$ performed up to $Q_{1}=Q_{2}=1$ GeV yields only around $90\%,80\%$ and $70\%$ of the total result, with relative contributions above $95\%$ not reached up to $Q_{1}=Q_{2}=1.8,2.5,$ and $3$ GeV, respectively. The discussion above implies the following requirements for a precise calculation (e.g. below $10\%$) when describing the TFFs:

• •

An accurate, precise and ideally model-independent method which can be improved upon via including new theoretical constraints and new experimental data.

• •

The method should implement a full-energy TFF description for the whole SL region (the time-like (TL) region is not involved in Eq. (5)), including well-known low- and high-energy constraints, the former due to integral weights at low-energies, the latter to render the loop integrals finite.

• •

The method should provide a very precise description at energies as low as $1$ GeV and, at least, a precise description for higher energies up to around $2-3$ GeV.

We believe that none of the current approaches for describing the TFFs fulfill all the criteria enumerated above and an alternative approach is desirable if our goal is a $10\%$ error. The pioneering works in Refs. Bijnens et al. (1996, 2002); Hayakawa et al. (1995), which were based on large-$N_{c}$ or vector meson dominance (VMD) approaches Bijnens et al. (1996, 2002); Hayakawa et al. (1995), consisted on a model of the large-$N_{c}$ limit of QCD. As such, a typical large-$N_{c}$ error estimate was given to be $30\%$, which represented an adequate error, but it is not enough at the present requested precision. Their systematic uncertainties and achieved accuracy with respect to the real TFFs are difficult to ascribe or systematically improve. A possible venue to refine these approaches is to use them as fitting functions using the current large amount of existing SL data for the single-virtual TFFs (see for instance Refs. Knecht and Nyffeler (2002); Kampf and Novotny (2011); Roig et al. (2014) where such approach was pursued), which could endow them with certain model independency, or at least, an accuracy beyond the conventional $30\%$ estimate — a proof of concept is given in Appendix B. It is uncertain however up to which accuracy could these approaches describe the real TFFs since the models used to fit are valid only in the large-$N_{c}$ limit of QCD and, if precision requires, how to systematically improve them. Besides, and unfortunately, there is at present a lack of SL data below $0.8$ GeV 777The exception is the L3 data for the $\eta^{\prime}$ Acciarri et al. (1998) which, to our best knowledge, has never been used so far in $a_{\mu}^{\textrm{HLbL}}$ calculations., see Fig. 4 and Fig. 5.

As a result, their low-energy description does not rely on data fitting, but on a fit extrapolation. The precision that such extrapolation provides on the relevant low energies — even if they may provide an excellent description for the available SL data — is difficult to quantify. A possible estimate of the eventual precision reached can be obtained by comparing with the available low-energy TL data for the $\eta$ meson Arnaldi et al. (2009); Berghauser et al. (2011); Aguar-Bartolome et al. (2014); Arnaldi et al. (2016); Adlarson et al. (2016a); the accuracy achieved there should provide a reasonable estimate for their SL counterpart. The study performed in Ref. Escribano et al. (2015) suggests the presence of a non-negligible error for such extrapolations. The lack of ability of these approaches to precisely reproduce the single-virtual low-energy TL data could in addition suggest a similar or even larger uncertainty in the double-virtual region, where no data is available so far to constrain their reconstruction.

Summarizing, the present data suggest that the standard procedures and the reference studies are not optimal at low energies and, as we will justify in Sect. IV, they cannot be considered model independent and cannot be systematically improved upon up to an arbitrary precision.

More recently, a dispersive reconstruction for the pseudoscalar TFFs has been formulated in Refs. Hoferichter et al. (2014); Hanhart et al. (2013); Xiao et al. (2015). Such approach has the advantage of relying on a data-based framework. As a result, it could in principle be as precise as its required inputs are and its precision systematically improved accordingly. This apparently solves the weakness of previous approaches. However, its full-energy implementation is in practice involved and cannot be complete in a model-independent way. It is for this reason that it is in practice so far limited to the low-energy region up to around $1$ GeV Hoferichter et al. (2014) and does not incorporate the high-energy constraints in its present form — the difficulty only increases for the double-virtual description — the reasons for which we believe them insufficient to pin down the error on the SL integrals beyond $10\%$. Any improvement with respect to traditional approaches comes consequently at the cost of the necessary mid-, and high-energy descriptions. It would be desirable then an alternative approach able to deal with all experimental data and theoretical constraints in the SL region (find discussions in Appendix G).

It was proposed in Ref. Masjuan and Peris (2007) that VMD approaches or the minimal hadronic approximation could be understood, in the large-$N_{c}$ limit of QCD, within the mathematical framework of Padé approximants (PAs) for the case of meromorphic functions. Such framework would provide the desired systematic method to reconstruct the TFF up to arbitrary precision, improving in accuracy with respect to former VMD approaches. Not only this, PAs guarantee the appropriate low-energy behavior by construction, while allowing to implement, at the same time, the high energies. The application of PAs is however not restricted (as in VMD approaches) to the large-$N_{c}$ limit of QCD, but can be applied to real $(N_{c}=3)$ QCD quantities — a textbook example is the HVP Masjuan and Peris (2010); Aubin et al. (2012) — provided they meet certain analytic properties which, for instance, the dispersive approach in Ref. Hanhart et al. (2013) fulfills. As a consequence, they provide a complementary tool to dispersive approaches in the SL region for improving upon standard VMD descriptions at low energies, but with the advantage of retaining the appropriate mid- and high-energy behaviors as well. Finally, as a difference with respect to all the previous approaches, they provide a systematic error allowing for the desired model independency.

The PAs implementation for the single-virtual pseudoscalar TFFs was discussed for the first time, and in a data-driven way, in Refs. Masjuan (2012); Escribano et al. (2014, 2015, 2016). More recently, the excellent accuracy at low energies was proved when comparing to the recently released low-energy TL data for the $\eta$ and $\eta^{\prime}$ mesons Escribano et al. (2015, 2016), which corroborated the appropriate description at low energies. This implies in addition that, being based on analytic properties (the latter even allows to reproduce some results in which discontinuities are involved Masjuan and Sanchez-Puertas (2016)), our description would provide the appropriate extrapolation to the low-energy double-virtual region. However, describing the more general double-virtual TFF requires the extension of Padé approximants to the multivariate case and involves the use of Canterbury approximants (CAs), which main concepts are illustrated in the following section. This mathematical framework will be the basis for reconstructing the double-virtual TFFs required to calculate the $a_{\mu}^{\textrm{HLbL;}P}$ in a data-driven way in Sect. V.

IV Canterbury approximants

IV.1 Definitions

Given an analytic function, symmetric in its variables, $f(x,y)=f(y,x)$, and with a known formal series expansion,

[TABLE]

Canterbury Approximants (CAs) Chisholm (1973); Chisholm and McEwan (1974); Hughes Jones (1976) are defined as rational functions of polynomials $R_{N}(x,y)$ and $Q_{M}(x,y)$

[TABLE]

($b_{0,0}=1$ can be generally taken), which coefficients ($i\geq j$) $a_{i,j}\in\mathcal{N}$, and $b_{i,j}\in\mathcal{D}$ are defined as to satisfy the accuracy-through-order conditions Masjuan and Sanchez-Puertas (2016, 2015), i.e.,

[TABLE]

where $\textrm{dim}(\mathcal{E})=\textrm{dim}(\mathcal{N})+\textrm{dim}(\mathcal{D})$. This is, they are required to match certain terms of the (low-energy) series expansion in Eq. (12), which guarantees the correct and desired behavior at low energies. Besides, the high-energy constraints can be included by invoking a two-point CA, which is, as it is standard Baker and Graves-Morris (1996), requiring a set of accuracy-through-order conditions with the $f(x,y)$ expansion for $x,y\to\infty$888For a more detailed discussion and examples on how to reconstruct CAs and their performance, we refer the interested reader to the Appendix of Ref. Masjuan and Sanchez-Puertas (2016) — an additional exercise on convergence is provided in this section below.. The CAs approach guarantees, among others, the convergence in the cut complex plane of the $C^{N}_{N+1}(x,y)$ and $C^{N}_{N}(x,y)$ as $N\to\infty$999This means, everywhere except for the TL region above threshold production. Fortunately, this is irrelevant for $g_{\mu}-2$. Nevertheless, find comments when the TL region is involved in Refs. Escribano et al. (2015, 2016); Masjuan and Sanchez-Puertas (2016). sequences to meromorphic Cuyt (1990) and Stieltjes functions Alabiso and Butera (1975) (see Appendix F for definitions) , which are not only justified in the large-$N_{c}$ limit of QCD, but also in the light of the dispersive approach from Ref. Hanhart et al. (2013), respectively. Note however that the reconstruction outlined above forbids to identify the CA’s poles to the physical ones, and rules out typical VMD constructions as a systematic description for the pseudoscalar TFFs, even in the SL region, as previously anticipated101010For a reconstruction of this kind when the function is known to be meromorphic, we refer the interested reader to Ref. Masjuan and Peris (2007)..

IV.2 Toy Models

In order to show the expected performance of CAs, and to provide stronger confidence in our approach, we illustrate its operation, prior to the real case discussion, with the aid of two well-motivated but analytically very different models for the $\pi^{0}$ TFF111111Our interest in these models resides in the fact that they cannot be well described with a finite set of resonances and the power of PAs is highlighted. The complexity of QCD suggests as well these models to be realistic enough to capture the main QCD ingredients. On top, the convergence of the PA sequence is slow enough to be appreciated in numerical examples.. The first of them is a large-$N_{c}$ Regge model Ruiz Arriola and Broniowski (2006, 2010), and the second one is an extension of the logarithmic model from Ref. Radyushkin (2009) to the double-virtual case (see Appendix B for a detailed description). For both models convergence is expected since they belong to the class of meromorphic and Stieltjes functions, respectively; it is equally interesting however to test on the convergence rate, a property as relevant for us as convergence itself.

In the following, we choose the $C^{N}_{N+1}(Q_{1}^{2},Q_{2}^{2})$ sequence for evaluating $a_{\mu}^{\textrm{HLbL};\pi^{0}}$; the latter is not only motivated because of convergence properties, but also due to the loop integral in Eq. (5), which requires a $C^{N}_{M}(Q_{1}^{2},Q_{2}^{2})$ sequence for which $M>N$ if integration is to be taken up to infinity.

To illustrate the different possibilities that CAs offer, we proceed as follows. First of all, we reconstruct the TFFs models employing only the information which is provided from the low-energy TFF expansion, Eq. (12), alone. This reconstruction does not incorporate however the high-energy constraints which can be expressed as121212For the physical TFF, it is possible to implement the single-virtual high-energy asymptotic or Brodsky-Lepage (BL) behavior as well, see Sect. V. Unfortunately, this is not possible for the present models, which have a logarithmic enhancement — see Appendix B. In any case, we note that the $C^{N}_{N+1}(Q_{1}^{2},Q_{2}^{2})$ approximant already implements the correct $Q^{-2}$ BL behavior.

[TABLE]

Indeed, the chosen sequence behaves for large energies as $C^{N}_{N+1}(Q^{2},Q^{2})\sim Q^{-4}$ rather than $Q^{-2}$. Consequently, we implement in a second step the conditions above in a sequential manner, starting with the $Q^{-2}$ high-energy behavior (but not constraining the particular $\mathcal{C}_{1}$ value above) and progressively including additional $\mathcal{C}_{i}$ coefficients.

The results for the approximants, which are shown in Table 1 along with the models’ exact result, illustrate the following features

• •

The first approach (LE row of Table 1) shows a clear convergence according to our expectations; indeed, the second $N=1$ element already provides a small systematic error, well below $10\%$, for both models.

• •

Implementing the $Q^{-2}$ high-energy behavior (OPE0 row) shows a clear improvement on convergence, which is expected given the high-energy behavior of the integrand in Eq. (5).

• •

Further constraining the $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ coefficients (OPE1 and OPE2 rows) does not seem to affect or change convergence, which points to the relevance of the low energies and to the fast convergence of the method131313For the logarithmic model, $F_{P\gamma^{*}\gamma^{*}}^{\textrm{log}}(Q^{2},Q^{2})=F_{P\gamma\gamma}\frac{M^{2}}{M^{2}+Q^{2}}$. Consequently, all the OPE coefficients, $\mathcal{C}_{n}$, are trivially satisfied within our approach..

• •

Besides the particular systematic error for each model, we observe an expected more general feature: the systematic error of a given element $N$ can be inferred from its difference with respect to the $N-1$ element. This provides a model-independent estimation for the systematic uncertainty and, thereby, the sought model-independent result.

In addition, it is worth to comment on factorization approaches for which $C^{N}_{N+1}(Q_{1}^{2},Q_{2}^{2})\sim C^{N}_{N+1}(Q_{1}^{2},0)\times C^{N}_{N+1}(Q_{2}^{2},0)$. These are very popular and seem to represent a good approximation at low energies Maris and Tandy (2002); Bijnens et al. (2012); Masjuan and Sanchez-Puertas (2015); Xiao et al. (2015) (note that non-factorizable effects are formally of order $(Q^{2})^{2}$ in the low-energy expansion in any case). The results are shown in the sixth row of Table 1 (Fact row) and show a potential large systematic error. The latter is however not only due to the wrong behavior at high-energies — our low-energy reconstruction in the second row of Table 1 (LE row) does not fulfill it either — but to the fact that not even the series expansion factorizes.

Finally, in our discussion above, it cannot be overemphasized the relevance of having employed the low-energy expansion Eq. (12) when reconstructing the approximants — as the framework requires — rather than fitting the rational functions to data themselves. To illustrate this statement, we show in the last row of Table 1 (Fit${}^{\textrm{OPE}}$) what would have been obtained if fitting the $C^{N}_{N+1}(Q_{1}^{2},Q_{2}^{2})$ rational functions, with the OPE behavior implemented, to a double-virtual data grid ranging from $0\leq Q_{1,2}^{2}\leq 35~{}\textrm{GeV}^{2}$ with a $2.3~{}\textrm{GeV}^{2}$ spacing. The obtained convergence is slower, and illustrates the difference and the power of CAs with respect to standard fitting approaches.

Summarizing the previous results: we find that the $C^{N}_{N+1}(Q_{1}^{2},Q_{2}^{2})$ sequence of approximants provides an excellent convergence when calculating $a_{\mu}^{\textrm{HLbL};\pi^{0}}$ for the chosen TFF models — which is further accelerated if the high-energy behavior is accounted for. More important, we find that the systematic uncertainty can be estimated from the difference among the elements within the sequence, which represents the main advantage from our approach and provides for the model independency. Having introduced CAs, motivated a sequence and illustrated its performance, we proceed to apply this approach for the real QCD case.

V Results

For the physical TFF, we define the formal series expansion, Eq. (12), in terms of the low-energy parameters (LEPs) $b_{P},c_{P},a_{P;1,1},...$ as

[TABLE]

It turns out that, under certain approximations, the authors of Ref. Stollenwerk et al. (2012) proved the isovector contribution to the TFF to be a Stieltjes function (cf. Appendix F), for which convergence of Padé approximants is guaranteed in advance.

Actually, Padé theory not only provides a convergence theorem for a sequence of PAs to Stieltjes functions, i.e., $\lim_{N,M\to\infty}P^{N}_{M}(s)-f(s)=0$, but also its rate of convergence Baker and Graves-Morris (1996); Masjuan et al. (2008a); Masjuan and Peris (2010), which is given by the difference of two consecutive elements in the PA sequence Masjuan et al. (2008b); Masjuan (2012); Escribano et al. (2014, 2015).

Furthermore, in the large-$N_{c}$ limit of QCD, the TFF becomes a meromorphic function, for which convergence is guaranteed as well Masjuan and Peris (2007, 2008). The sum rule approach employed in Ref. Klopot et al. (2012) for describing the TFF is again of the Stieltjes kind. Moreover, our experience from analyses of the TFF from the SL data Masjuan (2012); Escribano et al. (2014) and the excellent predictions achieved in the low-energy TL region Escribano et al. (2015, 2016) suggests that convergence to the TFF is at work and that its relevant analytical properties are retained. We understand that all these features hold for the double-virtual case too.

The available analytical information on the TFF is scarce though; at low energies $F_{P\gamma\gamma}(0,0)$ is theoretically related in the chiral (and large-$N_{c}$ for the $\eta$ and $\eta^{\prime}$) limit to the Adler Adler (1969)-Bell-Jackiw Bell and Jackiw (1969) anomaly, and can be expressed as141414$\mathcal{Q}$ stands for the charge matrix, $F$ is the pion decay constant in the chiral limit, and $\lambda_{\pi,\eta,\eta^{\prime}}=\lambda^{3,8,0}$ in the chiral limit with $\lambda^{a}$ the Gell-Mann matrices and $\lambda^{0}=\sqrt{2/3}~{}\mathds{1}_{3\times 3}$.

[TABLE]

This expression is, strictly speaking, valid only at the LO in both the chiral and the large-$N_{c}$ limits of QCD. Corrections to it which involve, among others, the $\eta-\eta^{\prime}$ mixing at the given order Bickert et al. (2016) are calculated in terms of unknown low-energy constants Goity et al. (2002). For this reason, we use instead the experimental results for $P\to\gamma\gamma$ decays in order to avoid model dependencies, which relation to $F_{P\gamma\gamma}(0,0)$ follows from $|F_{P\gamma\gamma}(0,0)|=\sqrt{\frac{64\pi}{(4\pi\alpha)^{2}}\frac{\Gamma(P\to\gamma\gamma)}{m_{P}^{3}}}$.

At small but finite virtualities, there are no further available theoretical predictions, and higher LEPs in Eq. (16) are theoretically unknown. Still, some of these LEPs were extracted for the single-virtual case in a data-driven approach using PAs Masjuan (2012); Escribano et al. (2014, 2015, 2016) which, as said, have proven extremely accurate when confronting them against the low-energy TL data. A similar procedure for the most general double-virtual case would be possible once double-virtual experimental data becomes available. At high energies, the TFF behavior can be theoretically described within pQCD. For the case of a single-virtual photon the TFF is known to behave according to the Brodsky-Lepage (BL) Lepage and Brodsky (1980) asymptotic behavior

[TABLE]

where $\pi_{\infty}$ depends on the pion decay constant, whereas $\eta_{\infty}$ and $\eta^{\prime}_{\infty}$ depend on the mixing parameters and the singlet-axial current running effects Escribano et al. (2015, 2016). For the double-virtual case, the TFF behavior at high energies is obtained from the operator product expansion (OPE) Novikov et al. (1984); Jegerlehner and Nyffeler (2009)

[TABLE]

where the numerical values for the parameters introduced in Eqs. (16, 18, 19) can be found in Table 6. Remarkably, Eqs. (18) and (19) ensure the convergence of the integrands in Eq. (5) and suggest the use of the $C^{N}_{N+1}(Q_{1}^{2},Q_{2}^{2})$ sequence explored in the previous section. Given the current available information on the double-virtual TFF, only the first two elements can be reconstructed. They are expressed as

[TABLE]

The connection to the LEPs from Eq. (16) is already visible in Eq. (20); the relation of the $\alpha_{i,j}$ and $\beta_{k,l}$ parameters in Eq. (21) to the LEPs is involved enough as not to fit in a single line. First, the single-virtual parameters $F_{P\gamma\gamma}(0,0),\alpha_{1},\beta_{1}$ and $\beta_{2}$ must be reconstructed. $F_{P\gamma\gamma}(0,0)$ is related, as mentioned, to the $P\to\gamma\gamma$ decays and can be extracted from the experimental values in Ref. Olive et al. (2014); $\alpha_{1},\beta_{1}$ and $\beta_{2}$ are related to the linear, quadratic and cubic terms in the single-virtual low-energy expansion, Eq. (16). These three parameters have been extracted from a data-driven approach in Refs. Escribano et al. (2014, 2015, 2016) for the $\eta$ and $\eta^{\prime}$, where they have been referred to as $b_{P}m_{P}^{-2},c_{P}m_{P}^{-4}$ and $d_{P}m_{P}^{-6}$, respectively. Alternatively, the cubic term can be traded for the BL asymptotic behavior, which is extremely convenient for the $\pi^{0}$ given the precise theoretical prediction (which contrasts with the $\eta$ and $\eta^{\prime}$ cases, see Table 6). Consequently, for the $\pi^{0}$, we employ the linear and quadratic terms determined in Ref. Masjuan (2012)151515In the near future, the data which are being analyzed at BES III Adl (2014) Collaboration in the low-energy SL region will allow for an accurate extraction of $d_{\pi}$. together with the BL prediction, which implies $\lim_{Q^{2}\to\infty}Q^{2}C^{1}_{2}(Q^{2},0)=2F_{\pi}=0.1884(3)$ GeV. It remains to determine the double-virtual parameters $\alpha_{1,1},\beta_{1,1}$ and $\beta_{1,2}$. For the $\pi^{0}$ case, it is possible to relate two of them to the high-energy expansion $P_{\infty}$ and $\delta_{P}$ parameters in Eq. (19). For the $\eta$ and $\eta^{\prime}$, $\delta_{P}$ is unknown; we take $\delta_{\eta,\eta^{\prime}}=\delta_{\pi}$ and ascribe an extra $30\%$ systematic error from $SU(3)_{F}$-breaking (and large-$N_{c}$) effects161616We note that such error covers for the observed $\pi^{0},\eta$ and $\eta^{\prime}$ differences for all the parameters which have been determined so far: $F_{P\gamma\gamma}(0,0),b_{P}m_{P}^{-2},P_{\infty},...$ . Besides, we find that, in practice, the $a_{\mu}^{\textrm{HLbL};P}$ dependence on this parameter is certainly mild.. Finally, one parameter remains to be determined. This can and should be related to the low-energy parameter $a_{P;1,1}$ in Eq. (16), which could be determined if double-virtual data becomes available. Hopefully, this may be possible in the future at BES III Adl (2014) for the $\pi^{0}$. Additional sources of information would be $P\to\bar{\ell}\ell$ Masjuan and Sanchez-Puertas (2015, 2016) and $P\to\bar{\ell}\ell\bar{\ell^{\prime}}\ell^{\prime}$ decays Escribano and Gonzalez-Solis (2015). Unfortunately, the available experimental precision for these measured decays is still not enough to make an extraction — find more comments later Section VI. For this reason, and to be as model independent as possible, we take for $a_{P;1,1}$ the most general range which is physically accessible without spoiling the TFF’s analytic properties, i.e., avoiding the presence of poles for the $C^{1}_{2}(Q_{1}^{2},Q^{2}_{2})$ in the SL region. This leads in practice to a range of the kind $a_{P;1,1}^{\textrm{min}}\leq a_{P;1,1}\leq a_{P;1,1}^{\textrm{max}}$171717Particularly, we find that $a_{\pi;1,1}^{\textrm{min}}=1.89b_{\pi}^{2}$, $a_{\eta;1,1}^{\textrm{min}}=1.65b_{\eta}^{2}$, and $a_{\eta^{\prime};1,1}^{\textrm{min}}=1.32b_{\eta^{\prime}}^{2}$, whereas $a_{\pi;1,1}^{\textrm{max}}=2.10b_{\pi}^{2}$, $a_{\eta;1,1}^{\textrm{max}}=6.00b_{\eta}^{2}$ and $a_{\eta^{\prime};1,1}^{\textrm{max}}=3.41b_{\eta^{\prime}}^{2}$. and completes the discussion about the CAs reconstruction.

The numerical integrals have been calculated with Mathematica 8.0 using the AdaptiveQuasiMonteCarlo method. The obtained results for the $C^{0}_{1}(Q_{1}^{2},Q_{2}^{2})$ are collected, for the different pseudoscalars, in Table 2 first column; the results for the $C^{1}_{2}(Q_{1}^{2},Q_{2}^{2})$ are collected, for the $a_{P;1,1}$ considered range, in the second and third columns of Table 2.

The errors are statistical only and arise from a MC analysis. For the $C^{0}_{1}$ case, the first and second errors arise from the TFF normalization ($F_{P\gamma\gamma}(0,0)$) and the slope parameter ($b_{P}$), respectively; for the $C^{1}_{2}$ $a_{P;1,1}^{\textrm{min}}$ choice (second column in Table 2), the error is due to the single virtual parameters $F_{P\gamma\gamma}(0,0),b_{P},c_{P},d_{P}$ and $P_{\infty}$181818In this limiting case, the $\alpha_{1,1}$ and $\beta_{2,1}$ parameters in Eq. (21) vanish, and the $\delta_{P}$ parameter becomes as a consequence irrelevant.; for the $C^{1}_{2}$ $a_{P;1,1}^{\textrm{max}}$ choice, an additional error related to $\delta_{P}$ arises; the last error in brackets stands in every case for the combination in quadrature of the individual ones.

The last line, Total, is the sum of the $\pi^{0},\eta$ and $\eta^{\prime}$ contributions. To estimate the systematic uncertainty, we take the largest difference among the $a_{P;1,1}^{\textrm{min}}$ and $a_{P;1,1}^{\textrm{max}}$ choices with respect to the $C^{0}_{1}(Q_{1}^{2},Q_{2}^{2})$, which even if it may overestimate the systematic uncertainty, still allows to reach an error below $10\%$. We note that this procedure corresponds to assume a fully correlated systematic error among the different pseudoscalars systematic uncertainties. Our final result is:

[TABLE]

where in the last line, the mean value among the third and fourth column in Table 2 has been taken, and their difference associated as an additional source of uncertainty which has been added in quadrature to the previous ones. In Appendix E we combine our final result with the Glasgow consensus in an attempt to provide a value for the whole $a_{\mu}^{\textrm{HLbL}}$.

Our obtained result can be compared to former determinations of the pseudoscalar-pole contribution, which appear in Refs. Knecht and Nyffeler (2002); Roig et al. (2014). The first of them Knecht and Nyffeler (2002), which was intended to clarify a sign discussion, reads $a_{\mu}^{\textrm{HLbL;}P}=(58(10)+13(1)+12(1)=83(12))\times 10^{-11}$ and has a crude error estimation as discussed in Ref. Knecht and Nyffeler (2002). Particularly, the $\eta$ and $\eta^{\prime}$ contributions were provided as an order of magnitude estimate, and for that reason only a simplified (and factorized) model was used. In addition, we note that our approach incorporates far more data for the TFFs which have become available since this study appeared.

The second and more recent study, Ref. Roig et al. (2014), is based on resonance chiral theory and obtains $a_{\mu}^{\textrm{HLbL;}P}=(57.5(6)+14.4(2.6)+10.8(0.9)=82.7(6.6))\times 10^{-11}$. Their analysis makes use of a similar data set to that employed in our approach for the $\pi^{0}$ case, but only up to 2014. The difference for the $\pi^{0}$ contribution could be ascribed to their different high-energy behavior since they cannot incorporate both BL and OPE, Eqs. (18), (19) — a known issue when using a single resonance Masjuan and Peris (2007). For the $\eta$ and $\eta^{\prime}$ cases, their approach description becomes more involved, as an appropriate $\eta-\eta^{\prime}$ mixing description requires an analysis at the next-to-leading order within their approach, with an error difficult to quantify. For this reason, Ref. Roig et al. (2014) uses instead $U(3)_{F}$ symmetry arguments to relate both $\eta$ and $\eta^{\prime}$ TFFs with the $\pi^{0}$ one. The differences we find illustrate the importance of a data-based approach to describe these TFFs.

As a final remark, both approaches in Refs. Knecht and Nyffeler (2002); Roig et al. (2014), which rely on the large-$N_{c}$ limit, can be understood as a CA to meromorphic functions. Then, the (missing) systematic error for their reconstruction — in which the poles are fixed in advance — is larger than in our case, see Ref. Masjuan and Peris (2007) and comments in Appendix B. Moreover, they do not employ a low-energy reconstruction, but a fitting procedure which, as illustrated in Sec. IV, entails an even larger error. Overall, these considerations suggest total unaccounted errors above $10\%$ for these approaches. An additional calculation based on the VMD models from Ref. Knecht and Nyffeler (2002) but fitted to lattice simulations for the $\pi^{0}$ TFF obtained $a_{\mu}^{\textrm{HLbL};\pi}=65.0(8.3)\times 10^{-11}$ Gérardin et al. (2016). The latter includes statistical error, but lacks the inclusion of a systematic uncertainty inherent to the large-$N_{c}$-based fitted model.

VI The role of future data

Finally, we outline the impact that future data would have in our TFFs reconstruction and thereby in our $a_{\mu}^{\textrm{HLbL;}P}$ determination. First of all, the most relevant parameters are the values of the TFFs at zero virtualities. To see this, note from Eq. (5) that the whole contribution is proportional to the square of this quantity. As such, a relative $\Delta_{F_{P\gamma\gamma}}$ error on the former directly translates into a $2\Delta_{F_{P\gamma\gamma}}$ relative error for $a_{\mu}^{\textrm{HLbL;}P}$. In this respect, and in the light of Table 2, a further reduction on this quantity will substantially improve our total error. The PrimEx-II experiment at JLab Gasparian (2016) and the experiment planned at KLOE-2 Babusci et al. (2012) on the $\gamma\gamma\rightarrow\pi^{0}\rightarrow\gamma\gamma$ reaction will halve the $\pi^{0}$ contribution error. Regarding the $\eta$ and $\eta^{\prime}$, the future GlueX experiment at JLab Dudek et al. (2012) would allow to reduce the $\eta$ and $\eta^{\prime}$ counterpart associated error — only a $3\%$ final precision for the $\eta$ case has been reported so far Gan (2015), which would again halve the current error.

Second, from Section V and the results in Table 2, we find pressing to get new data on the $\pi^{0}$ TFF. This would allow to improve our $b_{\pi}$ and $c_{\pi}$ determinations — especially the dominant systematic error, see Ref. Masjuan (2012) — and to extract the $d_{\pi}$ parameter, providing therefore an alternative single-virtual description in terms of LEPs alone and no high energy coefficients, which would be a valuable cross-check of our results. This will be possible in the near future once the ongoing analysis at BES III Adl (2014) Collaboration becomes published — in addition, further data is expected at low energies from KLOE-2 Babusci et al. (2012) and GlueX experiment Gan (2015).

In addition, in view of the discussion raised by BABAR Aubert et al. (2009) data concerning the TFF behavior at high energies, an interesting test could be to reconstruct the $\pi^{0}$ TFF from the light-quark content of the $\pi^{0}$ meson. From our knowledge of the $\eta$ and $\eta^{\prime}$ TFFs and the $\eta-\eta^{\prime}$ mixing in the flavor basis it is possible to extract from experimental data a pure light-quark TFF and use it to calculate what we denominate the light-quark $a_{\mu}^{\textrm{HLbL;}l.q.}$ which should be similar to the actual $\pi^{0}$ one. Appendix D contains the detail of such calculation and the results reported there support our final result from Eq. (22).

Furthermore, there are efforts to measure the $\eta$ and $\eta^{\prime}$ TFFs at low SL energies at BES III Adl (2014) — a similar study would be possible at the GlueX experiment as well Dudek et al. (2012); Gan (2015).

Finally, it remains an important task to get the first information on the double-virtual TFF given the $a_{P;1,1}$ parameter-induced error. This could be possible for the $\pi^{0}$ at BES III Adl (2014). Such a measurement would not only allow to improve on our current estimate, but to eventually obtain further double-virtual parameters. In addition, the latter would allow to trade the OPE expansion parameters in Eq. (19) in favor of the LEPs in Eq. (16) when reconstructing the TFF. As said, it is possible to use as well $P\to\bar{\ell}\ell$ decays. However, the current precision is not sufficient to provide a competitive constraint, see Appendix C for details. Similarly, $P\to\bar{\ell}\ell\bar{\ell^{\prime}}\ell^{\prime}$ decays provide an interesting potential source of information. Unfortunately, once more, a high precision is required for them in order to provide competitive constraints — around a $5\%$ precision for the yet not measured $\eta^{(\prime)}\to 2\mu^{+}2\mu^{-}$ decays, see Ref. Escribano and Gonzalez-Solis (2015). Alternatively, there is the possibility that, in the future, lattice simulations at the physical $\pi^{0}$ mass could provide valuable information in determining the LEPs in our method — see advances in Refs. Feng et al. (2012); Cohen et al. (2008); Lin and Cohen (2013); Gérardin et al. (2016).

Finally, our approach could incorporate the low-energy parameters predicted from dispersive theory Hanhart et al. (2013); Hoferichter et al. (2014); Xiao et al. (2015). Not only that, but our approach would allow to extend these frameworks at higher energies, providing an analytic continuation to the high-energy region, where these approaches cannot apply (find further comments in Appendix G).

VII Conclusions

In this work, we have employed the mathematical framework of Canterbury approximants (bivariate Padé approximants) in order to reconstruct the double-virtual pseudoscalar TFFs and calculate the pseudoscalar-pole contribution to $a_{\mu}^{\textrm{HLbL}}$. The method allows to incorporate, at the same time, both the low- and the high-energy information on the TFFs. Though the former play the most relevant role in this calculation, a precise calculation should consider both simultaneously. This feature represents the first advantage from our method with respect to resonance approaches or dispersive representations.

The required information for the reconstruction of CAs has been obtained from our works in Refs. Masjuan (2012); Escribano et al. (2014, 2015, 2016) and employs data from over 13 different collaborations. As a novelty of our approach, the method provides a systematic treatment — which advantage with respect to resonance approaches is especially obvious when dealing with the $\eta$ and $\eta^{\prime}$ — and allows for a systematic error estimation, which provides for the model independence of the result and the second advantage with respect to existing approaches.

As a result, we have found $a_{\mu}^{\textrm{HLbL};P}=94.3(5.3)\times 10^{-11}$, which is larger than previous estimates by a quantity which essentially corresponds to future experiments’ uncertainty. We note that such quantity is extremely interesting, as it is not only the dominating contribution to $a_{\mu}^{\textrm{HLbL}}$, but it is present in any data-driven approach for calculating $a_{\mu}^{\textrm{HLbL}}$ so far. Furthermore, our approach will benefit in the near future from the large amount of data which is expected to appear, including the SL one which, so far, is not included in dispersive approaches. Concerning the latter, our approach can also benefit from the TFF’s dispersive representations, as soon as the LEPs are reported, to account for the $\pi^{0},\eta$ and $\eta^{\prime}$ contributions to $a_{\mu}$, which shows the flexibility and complementarity of our approach

Acknowledgements.

The authors would like to thank A. Nyffeler and K. Kampf for discussions and to S. Leupold for comments concerning the dispersive representation. Work supported by the Czech Science Foundation (grant no. GACR 15-18080S).

Appendix A Functions involved in the HLbL

The pseudoscalar-pole contribution to the HLbL tensor is given as Knecht and Nyffeler (2002)

[TABLE]

where $\epsilon^{0123}=+1$, $q_{1},q_{2},q_{3}$ are outgoing from the blob depicted in Fig. 1 and $k=q_{1}+q_{2}+q_{3}$ is incoming to it.

The $I_{1,2}(Q_{1},Q_{2},t)$ functions involved in the $a_{\mu}^{\textrm{HLbL};P}$ calculation, Eq. (5), are defined as

[TABLE]

[TABLE]

where the following functions have been employed

[TABLE]

Appendix B TFF models

B.1 Regge model

In this appendix we discuss in detail the Regge model introduced in section IV and compare it with an approximation to it based on a finite set of resonances. Our goal is to show how an approximation built from a finite set of resonances converges to the original model, which contains on turn an infinite set of them. Several strategies are explored. The Regge model in question is taken from Refs. Ruiz Arriola and Broniowski (2006, 2010) and reads

[TABLE]

where $\psi^{(n)}(z)=\partial^{n}\ln\Gamma(z)$ is the polygamma function, the parameter $a$ is the string tension (which is fixed to $1.3$ GeV2 based on the study from Ref. Masjuan et al. (2012)) and $M=0.708$ GeV is chosen to match the slope parameter $b_{\pi}$ Masjuan (2012). At high energies it behaves as $Q^{-2}\ln Q^{2}$ for the single virtual case Ruiz Arriola and Broniowski (2006), whereas it behaves as $aF_{\pi^{0}\gamma\gamma}/(\psi^{(1)}(M^{2}/a))Q^{-2}$ for $Q_{1}^{2}=Q_{2}^{2}\equiv Q^{2}$ when $Q^{2}\to\infty$.

The Regge model can be expressed as well as as an infinite sum over the resonances within the Regge trajectory weighted by their correspondent residues,

[TABLE]

A resonance approach to this model will consist in retaining a finite number of resonances, achieving an increased precision as soon as more terms in the sum are included, assuming of course that the model parameters (masses and residues) are known. Such an approximation has, however, a slow rate of convergence, which is well understood from Padé theory Masjuan and Peris (2007). To illustrate this, we perform a numerical test and show in the first row in Table 3 (called Res) what would be obtained for $a_{\mu}^{\textrm{HLbL};\pi^{0}}$ if truncating the sum in Eq. (30) for a finite number of resonances $n$ (e.g., up to $m=n-1$ in Eq. (30)).

The slow asymptotic convergence is attributed to the fact that not even the TFF at the origin, $Q_{1}^{2}=Q_{2}^{2}=0$, is precisely reproduced. Therefore, to improve on that, we do not use the residue of the heavier resonance in the truncation; instead, we choose to fix such parameter to reproduce the TFF at the origin. As shown in the third row of Table 3 (called Norm) this strategy improves considerably on convergence. Finally, we choose to match not one, but all the residues in the summation, to fulfill the low-energy expansion. The results are shown in the fourth row of Table 3 (called Der) and yield the expected improvement on convergence.

Besides, as customary in resonance approaches, one can fit the residues to a set of pseudodata instead of matching them to the low-energy expansion of the model. Performing such fit using the same points as for that in Table 1 Fit${}^{\textrm{OPE}}$ row, we obtain the results of fifth row in Table 3 (called Fit), which show an irregular convergence, if it converges at all. At this point, one could blame the incorrect behavior of the resulting approximant when both $Q_{1}^{2}=Q_{2}^{2}\equiv Q^{2}\to\infty$, that behaves as $Q^{-4}$ instead of $Q^{-2}$. Including additional terms to fulfill this behavior, we obtain the results in the sixth row in Table 3 (called Fit${}^{\textrm{OPE}}$), which show an improved convergence, but still not such a good convergence as CAs to the same model with the same pseudodata fit (cf. Table 1). This appendix illustrates the potential systematic errors when the poles in rational approximants are fixed in advance.

B.2 Logarithmic model

Besides the Regge model, we introduce the logarithmic one employed as well in Section IV. The latter is inspired in flat distribution amplitudes models as introduced in Ref. Radyushkin (2009), and is extended to the double-virtual case as follows

[TABLE]

where $M=0.530$ GeV is chosen again to reproduce $b_{\pi}$ Masjuan et al. (2012). Again, it is straightforward to see its large $Q^{2}$-behavior $Q^{-2}\ln Q^{2}$ for a single-virtual photon; for equal virtualities, $F_{P\gamma^{*}\gamma^{*}}^{\textrm{log}}(Q^{2},Q^{2})=F_{P\gamma\gamma}M^{2}(M^{2}+Q^{2})^{-1}$ and, as a consequence $\lim_{Q^{2}\to\infty}F_{P\gamma^{*}\gamma^{*}}^{\textrm{log}}(Q^{2},Q^{2})=F_{P\gamma\gamma}M^{2}Q^{-2}$.

Appendix C $P\to\bar{\ell}\ell$ decays impact

A further possibility to test the TFF double-virtual behavior is given by the $P\to\bar{\ell}\ell$ decays. Whereas these decays could never be directly employed to extract the TFF $q^{2}$-dependence, as it happens with $P\to\bar{\ell}\ell\bar{\ell}^{\prime}\ell^{\prime}$ decays, they offer an indirect probe in terms of a loop integral over the double-virtual TFF Masjuan and Sanchez-Puertas (2015). The proposal to use these decays as a constraint for $a_{\mu}^{\textrm{HLbL;}P}$ was considered for the first time in Ref. Hayakawa and Kinoshita (1998), but, to our best knowledge, it has not been seriously considered so far.

In Refs. Masjuan and Sanchez-Puertas (2015, 2016), we performed a detailed and careful study of these decays employing our method of CAs. So far, only the $\pi^{0}\to e^{+}e^{-}$ Abouzaid et al. (2007) and $\eta\to\mu^{+}\mu^{-}$ Abegg et al. (1994) have been measured. For the $\pi^{0}$, we found its experimental value $2\sigma$ away from our CAs prediction Masjuan and Sanchez-Puertas (2015). Particularly, we found that in order to reproduce the experimental result would require $\delta_{\pi}^{2}\gtrsim 10~{}\textrm{GeV}^{2}$ in Eq. (19) and $a_{\pi;1,1}<-4b_{\pi}^{2}$ in Eq. (16). This would imply large corrections to the leading OPE behavior, a result far from theoretical expectations. In any case, taking these values for reconstructing the TFF191919Here we assume the absence of new-physics effects., we would obtain that $a_{\mu}^{\textrm{HLbL;}\pi}$ would shift down to around $36(7)\times 10^{-11}$ Masjuan and Sanchez-Puertas (2015). Of course, this result calls for a new experimental determination of the $\pi^{0}\to e^{+}e^{-}$ branching ratio, which could be possible at the NA62 experiment. Note that a first measurement of the double-virtual TFF at BES III would discard the $\pi^{0}$ TFF as the explanation for the measured branching ratio. Looking forward into the future, we study here the precision that such experiment would require in order to provide a valuable constraint to our study; this is, to improve our current $1.89b_{\pi}^{2}\leq a_{\pi;1,1}\leq 2.10b_{\pi}^{2}$ range. We find that, due to the narrow range for $a_{\pi;1,1}$, a precision of $0.1\%$ on the BR would be required. Still, even if this precision may be out of experimental reach, a new experiment that would shed light on the nature of the current deviation is highly desirable.

For the $\eta$ case, the sensitivity to the $a_{\eta;1,1}$ LEP is even lower. Even if an accurate description requires a proper double-virtual implementation, after implementing the parameters in Eq. (19), there is not much sensitivity to the $a_{\eta;1,1}$ parameter. Particularly, a precision below $0.1\%$ in the branching ratio measurement would be required to discern values within our range. Unfortunately, there are additional effects which are of further relevance at this precision, such as the implementation within our framework of the $\pi\pi$ cut, which would demand a more refined theoretical study Masjuan and Sanchez-Puertas (2016). In any case, a $2\%$ accuracy would be interesting already in order to corroborate (or falsify) our predictions, which may provide the only (indirect) experimental test of our TFF description in the near future. At present, the current precision is of $14\%$ and the central value is $1.4\sigma$ above our prediction. Similar comments apply for the $\eta^{\prime}$ decays too, which however have never been observed, being upper bounds the only available information so far Akhmetshin et al. (2015); Achasov et al. (2015). It would be therefore an interesting possibility in the future to access $\eta(\eta^{\prime})\to\mu^{+}\mu^{-}$ decays at LHCb Huong et al. (2016), as well as the possibility to measure the $e^{+}e^{-}\to\eta^{\prime}$ process at KLOE-II Solodov .

Appendix D The light-quark TFF

It is widely assumed that the smallness of the OZI rule violation would allow to express the $\eta$ and $\eta^{\prime}$ TFFs in terms of the light- and strange-quark TFFs Feldmann (2000). In addition, it is believed that the former would resemble the $\pi^{0}$ one except for a $5/3$ charge factor. Under the validity of these assumptions, the knowledge of the $\eta-\eta^{\prime}$ mixing provides thereby a cross-check for the $\pi^{0}$ TFF and, consequently, $a_{\mu}^{\textrm{HLbL};\pi^{0}}$.

In our analyses from Refs. Escribano et al. (2014, 2015, 2016), we could determine the mixing parameters, obtaining small violations of the OZI rule and confirming thereby the accuracy of such an approximation202020From these studies we obtain the mixing angles $\phi_{q}\simeq\phi_{s}\simeq 39.5^{\circ}$.. From these results and those for $b_{\eta^{(\prime)}},c_{\eta^{(\prime)}},d_{\eta^{(\prime)}}$ and $\eta^{(\prime)}_{\infty}$ (cf. Table 6), we obtain a light-quark $\eta_{q}$ TFF which resembles the $\pi^{0}$ in the SL region, and seems to support Belle tendency Uehara et al. (2012) against the BABAR Aubert et al. (2009) one at high energies, see Fig. 6. For completeness we also show the strange-quark TFF understood as the one for a pure $s\bar{s}$ $\eta_{s}$ state. Fig 6 also shows the experimental data for the $\pi^{0}$ TFF from Refs. Olive et al. (2014); Behrend et al. (1991); Gronberg et al. (1998); Aubert et al. (2009); Uehara et al. (2012); del Amo Sanchez et al. (2011); Acciarri et al. (1998); del Amo Sanchez et al. (2011) (blue circles) and the experimental data for the $\eta$ and $\eta^{\prime}$ TFFs rotated to the flavor light-strange-quark basis. The resulting parameters from the $\eta_{q}$ TFF read

[TABLE]

Indeed, every parameter is compatible with those of the $\pi^{0}$ TFF except for the normalization.

Taking the assumption that $F_{\pi^{0}\gamma^{*}\gamma^{*}}(Q_{1}^{2},Q_{2}^{2})=(3/5)F_{\eta_{q}\gamma^{*}\gamma^{*}}(Q_{1}^{2},Q_{2}^{2})$, we can provide an alternative determination for $a_{\mu}^{\textrm{HLbL;}\pi}$. This is found in the second row of Table 4, called l.q. A non-negligible shift appears with respect to our $\pi^{0}$ results in Table 2. However, this is to be expected given their different normalizations, which in turn represent the most relevant parameter in the calculation. Normalizing the light-quark TFF to match the $\pi^{0}$ one, an excellent agreement is found as shown in the third row of Table 4, called l.q. norm.

The errors identification is analog to that in Table 2. The values obtained in this exercise reassess our main results in Table 2 and supports our statement that future experimental data at high SL energies would not change much the $a_{\mu}^{\textrm{HLbL};\pi^{0}}$ central value, whereas they will provide a preciser determination for the LEPs in Eq. (16).

Appendix E Beyond pole contribution

There is at present a debate on how to deal with the HLbL tensor high-energy behavior dictated by the OPE (cf. discussion in Jegerlehner and Nyffeler (2009) with respect to Melnikov and Vainshtein (2004) and the summary talk by Vainshtein in D’Ambrosio et al. (2016)). Whereas we do not want to enter this debate, in this appendix we discuss how our approach could be used by both approaches Melnikov and Vainshtein (2004); Jegerlehner and Nyffeler (2009). The first approach Melnikov and Vainshtein (2004) proposes to modify the $\pi^{0}$-pole contribution to $a_{\mu}$ as such that certain OPE constraint to the $\langle VVVV\rangle$ Green’s function is satisfied. Its modification results on setting the external TFF (the gray blob connected to the external photon in Fig. 2) to a constant one. Following such prescription and using our description for the pseudoscalar TFFs, we obtain the results for $a_{\mu}^{\textrm{HLbL;}P}$ shown in Table 5.

The larger errors obtained now are proportional to the larger central values with, essentially, the same proportionality than our main result. Accounting for the systematic error as we did in Section V, we would obtain

[TABLE]

to be compared with the result from Ref. Melnikov and Vainshtein (2004) $a_{\mu}^{\textrm{HLbL}}=(76.5+18+18)\times 10^{-11}\to 114\times 10^{-11}$. This comparison illustrates again the potential large systematic errors — beyond $10\%$ — typical of resonance models. Actually, the result from Melnikov and Vainshtein (2004) was used in the Glasgow consensus to obtain the reference value $10.5(2.6)$. If we would replace their pseudoscalar-pole contribution by our Eq. 32, the final result would be

[TABLE]

one sigma larger, and in better agreement with ballpark estimates.

The second approach Jegerlehner and Nyffeler (2009) provides instead a model for the pseudoscalar contribution (not only the lightest pseudoscalar poles), related to the $\langle VVP\rangle$ Green’s function. It would be possible within our approach to reconstruct an analogous Green’s function. In this scenario to satisfy all the constraints imposed by the OPE, one should start directly with the $C^{1}_{2}(Q_{1}^{2},Q_{2}^{2},(Q_{1}+Q_{2})^{2})$ since the $N=0$ is too limited. Then, we cannot provide a systematic error and check on convergence. Besides, further constraints on this $\langle VVP\rangle$ Green function would be desired. For these reasons, we decide not to give a value for such scenario. In any case, we remind that this approach, as well as the previous one, were inspired in the $\pi^{0}$ TFF model in Ref. Knecht and Nyffeler (2002) which, as said, entails non-accounted systematic errors.

Appendix F Stieltjes functions

A function is said to be of the Stieltjes kind if it admits an integral representation Baker and Graves-Morris (1996)

[TABLE]

where $\phi(u)$ is any bounded and nondecreasing function Baker and Graves-Morris (1996). To see that such is the case for the isovector contribution to the TFF in Refs. Stollenwerk et al. (2012); Hanhart et al. (2013); Kubis and Plenter (2015), let $R=4m_{\pi}^{2}$, and define $d\phi(u)={\rm const.}\times\frac{q^{2}}{\pi}\frac{\operatorname{Im}F(1/u)}{u}$; making the change of variables $u=1/s$, Eq. (34) returns the once-substracted dispersive representation of the isovector contribution discussed in Ref. Stollenwerk et al. (2012), and also exploited in Refs. Hanhart et al. (2013); Kubis and Plenter (2015), once $\operatorname{Im}F(s)=\sigma^{3}(s)P(s)|F_{V}(s)|^{2}$ is identified. Since $\sigma(s)=\sqrt{1-4m_{\pi}^{2}/s}$, $P(s)$ is a linear polynomial with positive slope and $F_{V}(s)$ the $\pi^{\pm}$ vector FF, then $\operatorname{Im}F(s)$ is a positive function, the requirement of $\phi(u)$ to be nondecreasing is fulfilled and the convergence of PAs to the TFF is guaranteed.212121If the function $f(z)$ is a Stieltjes function, its $n^{th}$-subtracted version is a Stieltjes function as well Baker and Graves-Morris (1996).

Appendix G Dispersion Relations

In this appendix we develop our statements concerning potential drawbacks of dispersive approaches for extending the TFF representation into the SL region beyond energies of the order of $1$ GeV. For this purpose, we employ a simplified approach inspired from Ref. Hanhart et al. (2013). Specifically, we take the definition in Eq. (17) of that reference for the once-subtracted dispersion relation for the $\eta$ TFF, while we adopt a simpler but reasonable description for the $\pi^{\pm}$ vector form factor based on Refs. Gómez Dumm et al. (2000); Gómez Dumm and Roig (2013).

The result obtained for the TL region (from $q^{2}=0$ up to $q^{2}=m_{\eta}^{2}$), accessible in the $\eta\to\gamma\bar{\ell}\ell$ Dalitz decays, is illustrated in the upper panel of Fig. 7, and shows a nice agreement with existing data even though the approximations performed;

this nice overlap contrasts with the situation in the SL region at higher energies, which is illustrated in the lower panel of Fig. 7 and shows a clear deterioration above $1$ GeV. Of course, we stress that this is an oversimplified dispersive model and avoids for instance heavier resonances or inelasticities as done in Refs. Kubis and Plenter (2015); Xiao et al. (2015), but we hope it is enough to illustrate the expected features for large $Q^{2}$ values from dispersive representations, which are tightly related to the subtraction procedure.

In this respect, it would be interesting to use the series expansion at $Q^{2}=0$ that dispersive approaches can provide — given their reliability at low energies — to supply further information for the rational approach adopted in this work. Such procedure would provide a possibility to extend dispersive approaches up to arbitrarily large $Q^{2}$-values in the SL or, at the very least, to implement their LEPs.

Appendix H Input parameters

In this appendix we quote the different inputs used in our calculation. They are collected in Table 6 together with their original reference(s) and the data involved in the analysis for such determination. In addition, we quote in the last column those experiments which could improve the present values for the employed parameters. These parameters together with the definitions in Sections III and IV and Appendix A should allow the reader to reproduce all the results presented in this work.

We note that we have taken the opportunity to incorporate the recently released data from Refs. Arnaldi et al. (2016); Adlarson et al. (2016a)222222The experimental data from Ref. Adlarson et al. (2016a) supersedes the one used Aguar-Bartolome et al. (2014) in our previous work Escribano and Gonzalez-Solis (2015), which we remove in consequence. and Refs. Adlarson et al. (2016b); Lazzeroni et al. (2016)232323Ref. Lazzeroni et al. (2016) does not include the data points. Therefore, to include their analysis we used two different strategies: the first consists in generating a single point using their fit result, whereas the second consists in including their slope as a fitting parameter and leads to analog results. We thank M. Koval for discussions. regarding the $\eta$ and $\pi^{0}$ TFFs. For the $\eta$ case, we obtain very similar results to those in Ref. Escribano and Gonzalez-Solis (2015), with the advantage of reaching the $P^{3}_{3}$ approximant and obtaining a reduced statistical uncertainty due to the improvement in the experimental precision. For the $\pi^{0}$ case, the changes in central values are due to the updated $\Gamma(\pi^{0}\to\gamma\gamma)$ value from Ref. Olive et al. (2014) with respect to that used in Ref. Masjuan (2012), the inclusion of systematic uncertainties for Behrend et al. (1991) and, mildly, the data from Refs. Adlarson et al. (2016b); Lazzeroni et al. (2016). The improvement on systematics is low as a consequence of the limited and low TL $q^{2}$-range and we adopt those in Masjuan (2012). As an example, for the slope we obtain $b_{\pi}=0.0336(29),0.0321(13)$ and $0.0315(15)$ for the $P^{N}_{1}$, $P^{N}_{N}$ and $P^{N}_{N}$ sequence with the Brodsky-Lepage asymptotic behavior built-in (we reach up to the $N=6,2$ and $3$, element respectively). As usual, we take the average as our final result Masjuan (2012); Escribano et al. (2014); Escribano and Gonzalez-Solis (2015); Escribano et al. (2016), which is quoted in Table 6 and includes the systematic error, which has been combined in quadrature with the statistical one.

In addition, we give the $Q^{2}F_{\pi^{0}\gamma^{*}\gamma}(Q^{2})$ parameterization obtained from our the highest element reached within the $P^{N}_{N}$ sequence with the BL behavior built-in ($N=3$),

[TABLE]

which coefficients are shown in Table 7 (only the central values are given).