Lattice calculation of the pion transition form factor \pi^{0}\to\gamma^{*}\gamma^{*}

# Lattice calculation of the pion transition form factor π0→γ∗γ∗

Antoine Gérardin PRISMA Cluster of Excellence and Institut für Kernphysik, Johannes Gutenberg-Universität Mainz, 55099 Mainz, Germany    Harvey B. Meyer PRISMA Cluster of Excellence and Institut für Kernphysik, Johannes Gutenberg-Universität Mainz, 55099 Mainz, Germany Helmholtz Institute Mainz, Johannes Gutenberg-Universität Mainz, 55099 Mainz, Germany    Andreas Nyffeler PRISMA Cluster of Excellence and Institut für Kernphysik, Johannes Gutenberg-Universität Mainz, 55099 Mainz, Germany
###### Abstract

We calculate the transition form factor in lattice QCD with two flavors of quarks. Our main motivation is to provide the input to calculate the -pole contribution to hadronic light-by-light scattering in the muon , . We therefore focus on the region where both photons are spacelike up to virtualities of about , which has so far not been experimentally accessible. Results are obtained in the continuum at the physical pion mass by a combined extrapolation. We reproduce the prediction of the chiral anomaly for real photons with an accuracy of about . We also compare to various recently proposed models and find reasonable agreement for the parameters of some of these models with their phenomenological values. Finally, we use the parametrization of our lattice data by these models to calculate .

preprint: MITP/16-079

## I Introduction

The anomalous magnetic moment of the muon provides one of the most precise tests of the Standard Model of particle physics Jegerlehner:2009ry (); Miller:2012opa (). It is known to comparable precision in experiment Bennett:2006fi () and theory but the results disagree by about standard deviations Agashe:2014kda () depending on the theoretical estimate. To interpret this tension as a sign of new physics, improving the accuracy is of primary importance. On the experimental side, new experiments at Fermilab and J-PARC are expected to reduce the error by a factor of four future_g-2_exp (). Therefore, a corresponding theoretical effort is necessary to fully benefit from the increased experimental precision. The theory error of is dominated by hadronic contributions: the hadronic vacuum polarization (HVP) and hadronic light-by-light scattering (HLbL). The first contribution can be related to the cross section hadrons using a dispersion relation such that the estimate can, in principle, be improved by accumulating more data. Also, in recent years, more and more precise lattice QCD calculations of the HVP have become available but are not yet competitive with the dispersive approach DellaMorte:2011aa (); Boyle:2011hu (); Burger:2013jya (); Chakraborty:2016mwy (). However, the HLbL contribution to the muon cannot fully be related to direct experimental information and current determinations usually rely on model assumptions where systematic errors are difficult to estimate Prades:2009tw (); Jegerlehner:2009ry (); Bijnens:2015jqa (). However, recently a dispersive approach was proposed HLbL_DR () which relates the, presumably, numerically dominant pseudoscalar-pole contribution, as depicted in Fig. 1, and the pion loop in HLbL with on-shell intermediate pseudoscalar states to measurable form factors and cross sections with off-shell photons: and . Furthermore, increasingly realistic lattice calculations of the HLbL contribution to the muon have been carried out recently Blum:2014oka (); Blum:2015gfa (). Also, the hadronic light-by-light scattering amplitude per se has been calculated on the lattice in Green:2015sra ().

Within the dispersive framework, the pseudoscalar-pole contribution requires as hadronic input the transition form factor describing the interaction of an on-shell pseudoscalar meson, , with two off-shell photons with virtualities and . The HLbL contribution is then obtained by integrating some weight functions times the product of a single-virtual and a double-virtual transition form factor for spacelike momenta Jegerlehner:2009ry (). For the pion, the weight functions turn out to be peaked at low momenta such that the main contribution to arises from photon virtualities below  KN_02 (); Nyffeler:2016gnb (), a kinematical range accessible on the lattice.

The single-virtual transition form factor for the pion in the spacelike region has been measured experimentally by several collaborations Behrend:1990sr (); Gronberg:1997fj (); Aubert:2009mc (); Uehara:2012ag () in a wide kinematic range, although only for . More precise data down to are expected soon from BESIII Denig:2014mma (). There are currently no data available for the double-virtual transition form factor, a first measurement is planned at BESIII BESIII_double_virtual (). The double-virtual form factor has also been addressed on the lattice in Lin:2013im (). Finally, the authors of Feng:2012ck () also considered the double-virtual form factor at a single lattice spacing but focused their study on the pion decay , i.e. they were interested mostly in the behavior of the form factor at very low momenta. Transition form factors of mesons were first addressed in the context of the in Dudek:2006ut (); Chen:2016yau ().

Here we compute the transition form factor on the lattice in the kinematical region relevant to hadronic light-by-light scattering in the . Several lattice spacings and pion masses are used to extrapolate our results to the physical point. Our calculation involves several technical improvements over previous calculations.

This paper is structured as follows. In Sec. II, we give the precise definition of the transition form factor, describe its phenomenology and theoretical constraints from QCD and introduce the models whose functional form we will use to parametrize our lattice data. In Sec. III, we describe the methodology of the lattice calculation, including the analytic continuation, the required Wick contractions and the kinematic setup that we choose. In Sec. IV, the lattice calculation itself is presented, with the final result for the transition form factor presented in Sec. IV.4. Sec. V compares our fits to the lattice data with the available experimental and theoretical information on the pion transition form factor and in Sec. V.2 the pion-pole contribution to HLbL in the muon is evaluated with the form factor determined on the lattice. The paper ends with a summary of what has been achieved and an outlook on possible future improvements. Several appendixes contain some derivations and further discussions of some technical aspects, as well as tables with detailed results of the fits.

## Ii The pion transition form factor

In Minkowski spacetime, the transition form factor describing the interaction between a neutral pion and two off-shell photons is defined via the following matrix element111Equivalently, the form factor is given by

 Mμν(p,q1)=i∫d4xeiq1x⟨Ω|T{Jμ(x)Jν(0)}|π0(p)⟩=ϵμναβqα1qβ2Fπ0γ∗γ∗(q21,q22), (1)

where and are the photon momenta, the on-shell pion momentum, , is the hadronic component of the electromagnetic current and where we use the relativistic normalization of states . We use the mostly minus metric, , the axial current is given by with a Pauli matrix, and the phase of the one-pion state is fixed by with . In the chiral limit and at low energy, the form factor is constrained by the Adler-Bell-Jackiw (ABJ) anomaly Adler:1969gk (); Bell:1969ts (). At the physical pion mass, there are corrections due to quark mass effects which can be captured to a large extent by replacing the pion decay constant in the chiral limit by the pion decay constant obtained from charged pion decay Agashe:2014kda (). This leads to the following theoretical normalization of the form factor:

 \vspace−0.03cmFπ0γ∗γ∗(0,0)=14π2Fπ. (2)

At leading order in QED, one gets for the decay rate

 Γ(π0→γγ)=πα2em3π4F2π0γ∗γ∗(0,0), (3)

where is the fine structure constant. Together with Eq. (2) this reproduces quite well the measured decay width  eV  Agashe:2014kda (). The PDG average is dominated by the PrimEx experiment Larin:2010kq () where a precision of 2.8% has already been achieved and a further reduction of the error by a factor of two is expected soon. For a detailed comparison of theory and experiment at the level of a few percent, higher order quark mass and radiative corrections need to be taken into account, using chiral perturbation theory () together with some form of resonance estimates of the relevant low-energy constants pi0_gamma_gamma_theory (); Moussallam:1994xp ().

On the other hand, at large Euclidean (spacelike) momentum, the single-virtual form factor has been computed in the framework of factorization in QCD (operator-product expansion (OPE) on the light cone) with a perturbatively calculable hard-scattering part and a nonperturbative pion distribution amplitude. At leading order in , one finds the Brodsky-Lepage behavior BL_3_papers ()

 Fπ0γ∗γ∗(−Q2,0)−−−−→Q2→∞2FπQ2. (4)

In this formula, the prefactor should be taken with caution since its value actually depends on the shape of the pion distribution amplitude used in the calculation, which is usually modeled. When we impose below the Brodsky-Lepage behavior according to Eq. (4), we will only demand a falloff of the form factor, without insisting that the prefactor be reproduced exactly. On the experimental side, the single-virtual form factor has been measured for spacelike momenta in the range by CELLO Behrend:1990sr () and for by CLEO Gronberg:1997fj (). Later BABAR Aubert:2009mc () and Belle Uehara:2012ag () obtained results at larger momentum transfers both in the range . However their results differ significantly at large momenta: the results of BABAR showed an unexpected slower falloff of the single-virtual form factor, while the Belle data are compatible with a Brodsky-Lepage behavior. In any case, however, the data suggest that the asymptotic behavior is approached only at a momentum transfer above , outside the kinematical range considered in this paper. An analysis by BESIII Denig:2014mma () should be released soon which will cover the low-momentum region more relevant for the muon .

Finally, the double-virtual form factor where both momenta become simultaneously large has been computed using the OPE at short distances. In the chiral limit the result reads Nesterenko:1982dn (); Novikov:1983jt ()

 Fπ0γ∗γ∗(−Q2,−Q2)−−−−→Q2→∞2Fπ3[1Q2−89δ2Q4+O(1Q6)], (5)

where order corrections are neglected and the quantity parametrizes the higher-twist matrix element in the OPE and was estimated in Ref. Novikov:1983jt () using QCD sum rules. In the double virtual case, no experimental data exist yet but some results from the BESIII experiment are expected in the coming years in the range BESIII_double_virtual (); Nyffeler:2016gnb (). Therefore, the dependence of the double-virtual form factor in the kinematical range of interest for the computation of the hadronic light-by-light contribution to the muon is still unknown and the available estimates all rely on phenomenological models Jegerlehner:2009ry (); Bijnens:2015jqa (). The model parameters are either fixed using theoretical and experimental constraints from various sources or by fitting the experimental data of the single-virtual form factor and then extrapolating to the double-virtual case, i.e. by assuming a factorization of the form factor . However, this method might be unreliable and a model-independent theoretical estimate of the transition form factor from lattice QCD is highly desirable. Another a priori model-independent approach is the use of a dispersion relation for the form factor Hoferichter:2012pm (); DR_pion_TFF (), which is based on general properties of analyticity and unitarity. For the practical implementation, however, some assumptions and approximations need to be made.

Different phenomenological models have been proposed in the literature to describe the form factor in the whole kinematical range, see Ref. FF_reviews () and references therein. The simplest model is the vector meson dominance (VMD) model, where the form factor is given by

 FVMDπ0γ∗γ∗(q21,q22)=αM4V(M2V−q21)(M2V−q22), (6)

where to reproduce the anomaly constraint (2) and with usually set to the meson mass. We will, however, treat and as free model parameters in our fits to the lattice data below. The VMD model is compatible with the Brodsky-Lepage behavior (4) in the single-virtual case. However, it behaves as when both photons carry large virtualities and falls off faster than the OPE prediction (5). The second model considered in this paper is the lowest meson dominance (LMD) model Moussallam:1994xp (); Knecht:1999gb (), within the large- approximation to QCD, which can be parametrized as

 FLMDπ0γ∗γ∗(q21,q22)=αM4V+β(q21+q22)(M2V−q21)(M2V−q22). (7)

Again, one can set to recover the anomaly constraint. The form factor behaves as in the double-virtual case and for reproduces the leading OPE prediction, which is imposed in the original LMD model by construction. On the other hand, the model does not reproduce the Brodsky-Lepage behavior for the single-virtual form factor (4) but tends to a constant at large Euclidean momentum for the off-shell photon. The original LMD model has no free parameters, but we will treat and as free parameters in our fits below.

Finally, in Ref. Knecht:2001xc () the LMD+V model has been proposed as a refinement of the LMD model where a second vector resonance () is considered, see Ref. Nyffeler:2016gnb () for a recent brief review of the model. The LMD+V model can simultaneously fulfill the Brodsky-Lepage and the leading OPE behavior. Using a slightly different parametrization from Ref. Knecht:2001xc (), it can be written as

 FLMD+Vπ0γ∗γ∗(q21,q22)=˜h0q21q22(q21+q22)+˜h1(q21+q22)2+˜h2q21q22+˜h5M2V1M2V2(q21+q22)+αM4V1M4V2(M2V1−q21)(M2V2−q21)(M2V1−q22)(M2V2−q22). (8)

We have the relation , and between the above parametrization and the original model parameters (defined in the chiral limit) and (the latter parameters include corrections proportional to powers of the pion mass). In the LMD+V model proposed in Ref. Knecht:2001xc () only the parameters (or ) are treated as free parameters while the masses and are set equal to the physical masses of the and mesons. Furthermore the anomaly constraint is imposed, , as is the Brodsky-Lepage behavior which leads to . The form factor also has by construction the correct leading OPE behavior in the double-virtual case when both photons carry large Euclidean momenta by setting . As pointed out in Ref. MV_04 (), the parameter can be fixed by comparing with the subleading term in the OPE in Eq. (5). Finally the parameter has been determined in Ref. Knecht:2001xc () by a fit to the CLEO data Gronberg:1997fj () for the single-virtual form factor. One then obtains the model parameters

 ˜h2 = 0.327\leavevmode\nobreak GeV3,[¯h2=−4(M2V1+M2V2)+(16/9)δ2=−10.63\leavevmode\nobreak GeV2], (9) ˜h5 = (−0.166±0.006)\leavevmode\nobreak GeV,[¯h5=(6.93±0.26)\leavevmode\nobreak GeV4]. (10)

Following Ref. Moussallam:1994xp (), information on can also be obtained from the decay (assuming octet symmetry) which leads to the less precise determination  Knecht:2001xc (). In our fits below, we will in principle treat the parameters and the masses and as free parameters. The additional factors in the term with in the numerator in Eq. (8) will lead to more stable fits later.

A summary of the different asymptotic limits for each model and from the theory is given in Table. 1.

## Iii Methodology

From this section on, we use Euclidean notation by default. In particular, time evolution is governed by rather than , and . However the four-vectors and are always understood to be Minkowskian, i.e. .

### iii.1 Extraction of the form factor

Using the method introduced in Ji:2001wha (); Ji:2001nf (), and first implemented on the lattice in Dudek:2006ut (), one can show that the matrix element of Eq. (1) can be written in Euclidean spacetime as Feng:2012ck ()

 Mμν=(in0)MEμν,MEμν≡−∫dτeω1τ∫d3ze−i→q1→z⟨0|T{Jμ(→z,τ)Jν(→0,0)}|π(p)⟩, (11)

where is a real free parameter such that and denotes the number of temporal indices carried by the two vector currents. To obtain this formula, it is important to assume that so that the integration contour does not encounter a singularity, where one of the photons can mix with an on-shell particle. Therefore, one is led to consider the following three-point correlation function on the lattice

 (12)

where

 τ=ti−tf (13)

is the time separation between the two vector currents and

 tπ=min(tf−t0,ti−t0) (14)

is the minimal time separation between the pion interpolating operator and the two vector currents. Inserting a complete set of eigenstates, we obtain the following asymptotic behavior

 τ>0: C(3)μν(τ,tπ)−−−−→tπ→∞−Zπa32Eπ∑→z⟨0|Jμ(→z,τ)Jν(→0,0)|π(p)⟩e−i→q1→ze−Eπtπ, (15) τ<0: C(3)μν(τ,tπ)−−−−→tπ→∞−Zπa32Eπ∑→z⟨0|Jν(→0,−τ)Jμ(→z,0)|π(p)⟩e−i→q1→ze−Eπtπ, (16)

where is the overlap factor of our interpolating operator with the pion state222With the choice of made below in Eq. (23), the overlap is given by the partially conserved axial current (PCAC) relation, , where is the average quark mass. and the factor in the denominator comes from the relativistic normalization of states. The large time behavior of the three-point correlation function (12) ensures that the pion is on shell and that the excited states contribution in the pseudoscalar channel is small. Finally, the overlap and the pion mass are extracted from the two-point correlation function

 C(2)(t)=a3∑→x⟨P(→x,t)P†(→0,0)⟩e−i→p→x−−−→t→∞|Zπ|22Eπ(e−Eπt+e−Eπ(T−t)), (17)

where is the temporal extent of the lattice. It is convenient to remove the explicit pion energy time dependence in the three-point correlation function and to define

 Aμν(τ)=limtπ→+∞C(3)μν(τ,tπ)eEπtπ. (18)

Then, from Eq. (11), can be obtained via

 MEμν=2EπZπ(∫0−∞dτeω1τAμν(τ)e−Eπτ+∫∞0dτeω1τAμν(τ))=2EπZπ∫∞−∞dτeω1τ˜Aμν(τ), (19) ˜Aμν(τ)=limtπ→+∞eEπ(tf−t0)C(3)μν(τ,tπ)={Aμν(τ)\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak τ>0Aμν(τ)e−Eπτ\leavevmode\nobreak \leavevmode\nobreak \leavevmode\nobreak τ<0. (20)

The integral (19) is convergent as long as333The bound applies in infinite volume. In finite volume, the threshold can be at a slightly different energy than . : the three-point correlation function falls off with a factor , with the energy of a vector state. We point out that it is , rather than which is most directly related to the matrix element of interest ; see Appendix A for more details.

### iii.2 Kinematic setups

On the lattice, the momentum of the pion is set explicitly through the pseudoscalar interpolating operator used in Eq. (12) and its energy is then imposed by the on-shell condition. We are also free to choose one vector current spatial momentum (e.g. ), being determined by the momentum conservation . Finally, in Eq. (19), we can vary continuously , with determined by the energy conservation . Therefore, the kinematical range accessible on the lattice can be parametrized by

 q21 =ω21−→q21, q22 =(Eπ−ω1)2−(→p−→q1)2. (21)

Choosing the pion reference frame where , both photons have back-to-back spatial momenta () and . The kinematic range corresponding to different choices of is plotted in Fig. 2 for two different lattice resolutions. As explained below, from a numerical point of view, different momenta can be obtained without any new inversion of the Dirac operator. Therefore, this setup is adapted to study the form factors at large and, for each ensemble, the three-point correlation function has been computed up to momenta . Using the Lorentz structure of the form factor (see Eqs. (1), (11) and (19)), with one or more temporal indices vanishes and the spatial components can be written in the form

 Akl(τ)=−iqklA(τ),qkl≡ϵklαβqα1qβ2=mπϵkliqi1, (22)

where is a scalar under the spatial rotation group. From we define in the same way. Averaging over equivalent momenta through the cubic group, the statistic can be significantly increased. The total number of equivalent contributions for each value of is summarized in Table 2.

### iii.3 Correlation functions

We use the following (anti-Hermitian) interpolating operator for the neutral pion ,

 P(x)=¯¯¯u(x)γ5u(x)−¯¯¯d(x)γ5d(x)=¯¯¯¯ψ(x)γ5τ3ψ(x). (23)

At the quark level, the three-point correlation function receives three contributions,

 C(3)μν(τ,tπ)=Cconnμν(τ,tπ)+Cdisc1μν(τ,tπ)+Cdisc2μν(τ,tπ). (24)

Let , , , and and are the electromagnetic charges. Only up and down quark contributions are considered in this paper. If one uses two “local” vector currents,

 Jlμ(x)=∑fQf ¯¯¯¯ψf(x)γμψf(x), (25)

the connected contribution to the three-point correlation function reads

 Cconnμν(τ,tπ) =a6∑→x,→y⟨Jν(→0,tf)Jμ(→y,ti)P†(→x,t0)⟩ e−i→q1→y ei→p→x −tr[τ3Q2] a6∑→x,→y⟨ ¯¯¯¯ψf(z)γνψf(z)¯¯¯¯ψf(y)γμψf(y)¯¯¯¯ψf(x)γ5ψf(x)⟩ e−i→q1→y ei→p→x =2tr[τ3Q2] ⟨a6∑→x,→yReTr[G(x,z)γνG(z,y)γμG(y,x)γ5] e−i→q1→y ei→p→x ⟩U =2tr[τ3Q2] ⟨a3∑→yReTr[γνγ5G†(y,z)γ5γμ˜G(y,z;t0;→p)] e−i→q1→y⟩U, (26)

where Tr is the trace over spinor and color indices; is the Pauli matrix; is the charge matrix; tr the trace over flavors only; denotes the light quark propagator; and is a sequential propagator defined below. The correlation function depicted in Fig. 3 is computed in two steps: a first inversion on a point source leads to the solution vector . This solution vector is projected against the pion momenta and, when restricted to a given time slice , is used as a secondary source to obtain the sequential propagator (represented by a double line in Fig. 3)

 ˜G(y,z;t0;→p)=a3∑→xG(y,x)γ5G(x,z)ei→p→x. (27)

In particular, the sequential propagator satisfies the equation

 a∑yD(w,y)˜G(y,z;t0,→p)=δtw,t0γ5G(w,z)ei→p→w≡˜η(w), (28)

where is the sequential source and is the lattice Dirac operator. It is clear that a new sequential inversion would be required for each pion momentum and each value of . However, the momentum and the indices and can be chosen freely without any new inversion of the Dirac operator. This allows us to increase the statistics (see Table 2). In practice, we used ten sources per gauge configuration, randomly distributed on the lattice.

For the results presented in this paper, the connected three-point correlation function is computed using one local and one ‘point-split’ current. The latter is given by

 Jcμ(x) =∑fQf2(¯¯¯¯ψf(x+a^μ)(1+γμ)U†μ(x)ψf(x)−¯¯¯¯ψf(x)(1−γμ)Uμ(x)ψf(x+a^μ)). (29)

The Wick contraction is then only slightly modified. The point-split vector current satisfies the Ward identity and does not need any renormalization factor, contrary to the local vector current. In the -improved theory, the renormalized currents read

 Jα,Rμ(x)=ZαV(1+bαV(g0)amq)(Jαμ(x)+acαV∂νTμν), (30)

where the label stands for local or conserved and for isospin or , and are improvement coefficients and is the tensor density (written here for the improvement of the isovector part of the electromagnetic current). In particular, and , while the renormalization constant has been computed nonperturbatively in DellaMorteRD (); Fritzsch:2012wq () with a relative error below the percent level. In this paper we use the latter values both for the and currents. The improvement coefficients have been evaluated in Harris:2015vfa (), however in this study, we neglect the contribution from the tensor density as well as the improvement coefficient . Thus -improvement is only partially implemented.

For the disconnected contributions, we use two local vector currents. Wick contractions involving only the pion do not contribute since the and contributions exactly compensate each other. Therefore, one vector current must be contracted with the pion which leads to the two diagrams depicted in Fig. 4. The first diagram in the theory corresponds to the following contraction

 (τ,tπ)=−tr[Q] tr[τ3Q] a9V∑→x,→y,→z⟨ ¯¯¯¯ψf(z)γνψf(z)¯¯¯¯ψf(y)γμψf(y)¯¯¯¯ψf(x)γ5ψf(x)⟩U e−i→q2→z e−i→q1→y ei→p→x =−tr[Q] tr[τ3Q] ⟨a3∑→yTr[G(y,y)γμ]e−i→q1→y a6V∑→x,→zTr[G(z,x)γ5G(x,z)γν] e−i→q2→z ei→p→x ⟩U, (31)