Pion Transition Form Factor from Twisted-Mass Lattice QCD and the Hadronic Light-by-Light π 0 -pole Contribution to the Muon g − 2

The neutral pion generates the leading pole contribution to the hadronic light-by-light tensor, which is given in terms of the nonperturbative transition form factor F π 0 γγ ( q 21 , q 22 ). Here we present an ab-initio lattice calculation of this quantity in the continuum and at the physical point using twisted-mass lattice QCD. We report our results for the transition form factor parameterized using a model-independent conformal expansion valid for arbitrary space-like kinematics and compare it with experimental measurements of the single-virtual form factor, the two-photon decay width, and the slope parameter. We then use the transition form factors to compute the pion-pole contribution to the hadronic light-by-light scattering in the muon g − 2, finding a π 0 − pole µ = 56 . 7(3 . 2) × 10 − 11 .


I. INTRODUCTION
The anomalous magnetic moment of the muon (the muon g − 2) provides a stringent test of the Standard Model at high precision and the possibility to glimpse subtle effects of potential new physics beyond the Standard Model.Recent results from the Fermilab experiment [1,2] in combination with the Brookhaven result [3] have yielded the most precise experimental determination to date at the level of 19 ppm: a µ (exp) = 116 592 059(22) × 10 −11 , where a µ ≡ (g − 2) µ /2.A comparable precision is striven for on the theoretical side with the recent state of affairs summarized in the white paper [4].
The leading hadronic contributions to the muon anomalous magnetic moment come from diagrams involving the hadronic vacuum polarization (HVP) and the hadronic light-by-light (HLbL) scattering tensors.Both make significant contributions at the level of precision achieved by the experimental results.It is important to improve the theoretical determinations of both contributions to match future targets of experimental precision; this work addresses the HLbL contribution.
Two complementary approaches to the theoretical determination of the HLbL contribution are becoming recognizable for reliable estimates.On one hand, several direct lattice calculations have analyzed the complete HLbL tensor [5][6][7][8][9][10][11].On the other hand, the HLbL tensor can be decomposed into contributions from exchanges of various intermediate on-shell states in a data-driven dispersive approach, where contribu-P + flipped P FIG. 1. Pseudoscalar-pole diagrams contributing to the leading order hadronic light-by-light (HLbL )scattering.Each striped circle represents a nonperturbative P → γ * γ * transition form factor required to evaluate these contributions, where "P" indicates the exchanged pseudoscalar meson.The "flipped" contribution signifies the diagram with the pseudoscalar and photon portions of the diagram horizontally flipped.
tions from heavier intermediate states are suppressed [12][13][14][15].The latter approach has provided the most precise theoretical estimates of the HLbL contribution to the muon g − 2 to date [4], while the former has the benefit of allowing for a fully self-contained lattice QCD calculation from first principles without resorting to chiral effective field theory and experimental data.However, there has been an ongoing effort to provide also ab-initio data based on lattice QCD calculations for the leading contributions within the dispersive decomposition of the HLbL.
(3) Integrating over space-like q 1 , q 2 , and q 3 = −q 1 − q 2 with the appropriate kinematical factors from the diagrams in FIG. 1 yields contributions to the total muon g − 2 from each pseudoscalar meson exchange.
The π 0 -pole contribution to the muon g − 2 is larger than those of the η and η ′ by roughly a factor of four [17,18].This feature, together with the relative simplicity of accessing the pion TFF in both lattice and data-driven approaches, means that the study of this transition form factor is farthest advanced at this point.Direct experimental data for F π 0 γγ (q 2  1 , q 2 2 ) is primarily available in the singly virtual (either q 2 1 = 0 or q 2 2 = 0) regime with large photon virtuality [19][20][21][22][23]. Fits to this data using Canterbury approximants have yielded a data-based estimate across all virtualities [24].Meanwhile, the pion TFF has also been determined from indirect experimental data using a dispersive approach [25,26].On the other hand, an ab-initio determination is desirable to remove the dependence on experimental inputs.As such, the neutral pseudoscalar meson TFF has also been computed on the lattice.The calculation of the neutral-pion TFF on the lattice has been pioneered by the Mainz group using Wilson quarks [27,28].The BMW collaboration recently added their results using staggered quarks [29].Our preliminary results with Wilson twisted-mass quarks appeared in Refs.[30,31].The major sources of uncertainty in these approaches is the extrapolation to the continuum and, in case the quark masses are not tuned to their physical values as in Refs.[27,28], the extrapolation to the physical point.A first calculation of the η-meson TFF and the corresponding pole contribution at one lattice spacing has been reported by our collaboration in Ref. [32] and results for both the ηand η ′ -meson TFF and pole contributions addressing all systematic errors have been presented by the BMW collaboration in Ref. [29].
In the present work, we complement these existing efforts by evaluating the pion transition form factor in the continuum using the twisted-mass lattice QCD discretization, with ensembles generated directly at the physical point.This calculation serves to strengthen the lattice consensus of the continuum limit and minimizes systematic effects from extrapolation of the pion mass.In this work, we present a parameterization of the pion transition form factor and use it to calculate the π 0 -pole contribution to the muon g − 2. As our main result, we find a π 0 −pole µ = 56.7(3.2) × 10 −11 .We further use the transition form factor to estimate the two-photon decay width as Γ(π → γγ) = 7.50(0.50)eV, and the slope parameter as b π = 2.16(0.20)GeV −2 , which provides input for determining the electromagnetic interaction radius.Quoted uncertainties include statistical errors and systematic uncertainty from analysis choices and the continuum extrapolation.
The structure of this paper is as follows.In Sec.II, we review the theoretical background for the extraction of the required amplitude from lattice QCD simulations.In Sec.III, we describe the calculation of the transition form factors and their extrapolation to photon virtualities required for the observables that we calculate.In Sec.IV, we present the calculation of a π 0 −pole µ from the transition form factors and the determination of the decay width Γ(π → γγ) and the slope parameter b π .Next, in Sec.V, we discuss the estimation of statistical and systematic errors based on model averaging, before presenting our results in Sec.VI, comparing them to experimental and theoretical predictions.Then, in Sec.VIII, we make some concluding remarks and give an outlook on possible future directions for improving the calculation presented here.

A. Three-point function
We access the amplitude Ãµν (τ) through the following lattice three-point correlation function where P 0 † is an interpolating operator for the neutral pion and ⟨. . .⟩ indicates the Euclidean path-integral expectation value.Translation invariance, or equivalently conservation of momentum, ensures that the current j ν must carry momentum ⃗ q 2 = ⃗ p − ⃗ q 1 .We denote the minimal separation between the pion creation operator and the nearest electromagnetic current by t sep = min(τ, 0) + t π .When t sep ≫ 0, or equivalently t π ≫ 0, −τ, the pion creation operator P 0 † (⃗ p, −t π ) is in the distant Euclidean past, and the Euclidean matrix element in Eq. ( 4) can be recovered from the three-point function as where and ∆E indicates the energy gap to the first excited state produced by P 0 † (⃗ p, −t π ).Given this limiting behavior, we define In practice, we compute this quantity for finite values of t π taken sufficiently large to ensure the excited state contamination is suppressed below our uncertainties.When it is clear from context, we elide the t π argument in the following.
The parameters E π and Z π are extracted from a separate analysis of the the τ-dependence of the two-point function We note that for maximally twisted fermions the value of Z π is strictly positive and can be unambiguously extracted from Eq. ( 10), because it is related to the pion decay constant f π through Z π = f π m π sinh(m π )/2µ l where µ l is the twisted lightquark mass parameter [36,37].

B. Twisted-mass Lattice QCD
We perform our calculation using twisted-basis fermion fields, denoted χ.For the light quarks, these are related to the physical-basis doublet ψ = (u, d) through the maximal twist transformation ψ = e iπγ 5 τ 3 /4 χ , ψ = χ e iπγ 5 τ 3 /4 , (11) where τ 3 is the generator of the third component of isospin, represented by the third Pauli matrix.In the continuum limit, the twist transformation is a trivial change of basis.At finite lattice spacing, however, the Wilson term breaks chiral invariance meaning the twisted basis provides an alternative path towards the continuum limit, with certain desirable properties such as automatic O(a) improvement [36,38].
To relate calculations performed in the twisted basis to the physical amplitude and the TFF, we choose twisted basis operators which correspond to the physical-basis electromagnetic currents and pion operators.Specifically, we choose the pion annihilation operators to be where τ ± are the raising/lowering operators of isospin.We also choose local electromagnetic currents, given in the physical basis in the light-quark sector by where ) is the light-quark charge matrix.This current can be decomposed into terms of total isospin I = 1 and I = 0 (always with fixed z-component I z = 0), denoted respectively by j 1,0 and j 0,0 .The result is then easily written in the twisted basis, giving It will also be useful for the following section to define the currents with isospin I = 1 and I z = ±1, For the strange and charm quark we use a mixed-action approach: In the sea-quark action the strange and charm quark form a twisted-mass quark doublet, including tuned twistedmass and quark-mass splitting terms.Details about tuning of the twisted-mass parameters for the sea-quark action, in particular to maximal twist, are given in Ref. [39].In the valence action we use two Osterwalder-Seiler heavy-quark doublets [38] by introducing two doublets χ (h) = (χ (h)  + , χ (h) − ), for h = s, c, which avoids twisted-mass-related mixing of strange and charm flavor.Both doublets work analogously to the upand down-quark doublet, they have the same critical hopping parameter, but use individual tuning of the twisted quark-mass parameters.The strange-quark mass is tuned such that the Ω baryon takes its physical value, while the charm-quark mass is tuned to reproduce a near-physical Λ c mass.
The strange-and charm-flavor electromagnetic currents then use the local operators The iso-scalar strange-and charm-flavor currents enter correlation functions in FIG. 2 only as vector-current disconnected diagrams "V-disconnected I" and "V-disconnected II" (details on the Wick contractions are provided in the following section).To maintain automatic O(a) improvement we average the vector-current disconnected contributions from strange and charm propagators generated from both χ h + and χ h − with h = s, c.
We use gauge ensembles produced by the Extended Twisted Mass Collaboration (ETMC) from simulations of isospinsymmetric QCD (isoQCD) with N f = 2 + 1 + 1 flavors of Wilson Clover twisted mass quarks.The properties of the used ensembles most relevant for this work can be found in Tab.I; for more details see Refs.[39][40][41][42][43].
The local vector currents detailed above renormalize multiplicatively.Because of the twisted-mass formulation, isoscalar and neutral iso-vector vector currents are renormalized with renormalization factor Z V , while the charged iso-vector vector currents renormalize with Z A : These renormalization constants (RCs) were calculated to high precision for the ensembles used in this work in Ref. [43], as shown in Tab.II.The uncertainties on the RCs are negligible compared to the errors on the bare TFFs, such that we only use the central values of the RCs throughout this work.Two remarks are in order due to the fact that we are using nonconserved local vector currents for the construction of the three-point correlation function in Eq. ( 6).First, one can show that possible short-distance singularities are absent, i.e., the local currents in Eq. ( 6) admit a well defined continuum limit.The arguments are based on the operator product expansion and are given in Appendix D of Ref. [27] for Wilson fermions.Due to universality, they also apply to the Wilson twisted-mass lattice QCD formulation.Second, we note that the nonconserved local currents do not spoil the automatic O(a) improvement.This can be seen from the fact that all involved lattice quantities are constructed such that they have the correct symmetry property under the combined twisted-mass transformation P × D d × (µ → −µ) where P is the standard parity symmetry transformation, D d the operatordimensionality transformation, and µ the twisted-quark mass parameter.This symmetry forbids the appearance of O(a) terms in physical matrix elements for twisted-mass lattice QCD at maximal twist [36,38], and hence guarantees automatic O(a)-improvement of the three-point function in Eq. (6).
Connected II Wick contractions contributing to C ± µν (τ, t P ).These are the connected (top row), and vector-current disconnected ("Vdisconnected", bottom row) diagrams.Note that the connected diagrams only involve light quarks, while the disconnected vectorcurrent loop in the V-disconnected diagrams involves light, strange, and charm quarks for the 2 + 1 + 1 ensembles used in this work.

C. Wick contractions
A naive evaluation of the three-point amplitude Eq. ( 6) contains disconnected contributions in which the quarks within the neutral pion contract with themselves.These include a singly-disconnected contribution which correlates the neutralpion loop with the connected vector-vector two-point function, as well as a doubly-disconnected contribution with three loops.In the limit of exact isospin symmetry, these contributions vanish identically.In the twisted-mass formulation at finite lattice spacing, however, isospin symmetry is broken by O(a 2 ) lattice artifacts [37] (independently of the twist angle) and the disconnected contributions do not cancel at finite lattice spacing.
Since typically disconnected diagrams are costly to evaluate, we minimize the number of such diagrams by instead considering the isospin-rotated amplitude where spacetime coordinates have been suppressed for clarity.The amplitude in Eq. ( 20) does not contain pion-loop disconnected contributions, leaving only the connected and vector-current disconnected ("V-disconnected") contributions depicted in FIG. 2. This rotation is valid in the continuum due to the iso-symmetric continuum limit of twisted-mass lattice QCD.Consequently, the rotation only modifies the lattice artifacts at O(a 2 ) in the twisted-mass formulation, which is taken into account in the continuum extrapolation.We evaluate the connected diagrams using inversions based on point sources for the current j ν ( ⃗ 0, 0) with sequential inversion through the pseudoscalar operator P ± † .This means that each choice of momentum ⃗ p for the pseudoscalar requires a new inversion, while momentum ⃗ q 1 inserted at the sink current j µ (⃗ x, τ) can be easily varied.For the evaluation of the connected contribution to the three-point amplitude, 16 point sources per configuration are used for cB211.072.64 and 8 point sources per configuration on the other two ensembles cC211.06.80 and cD211.054.96 in Tab.I.
For the vector-current disconnected diagrams "Vdisconnected I/II" we use quark loops generated with hierarchical probing [44] in general, and low-mode deflation [45] in particular on ensemble cB211.072.64.Details on the dilution scheme for the Hadamard vectors and the number of low modes are given in Ref. [46] and in Table 2 of Ref. [47].The associated two-point correlators built from the charged pion and iso-vector vector current (cf.FIG. 2) are evaluated with stochastic timeslice propagators based on Z 2 × i Z 2 noise sources.The number of stochastic samples is chosen equal to the time extent in lattice units, T/a.

D. Kinematics
The pion momentum ⃗ p is set through the pseudoscalar interpolating operator while the momentum ⃗ q 1 is set through the current j µ , as shown in Eq. ( 6), with both momenta selected from the available finite-volume momenta The energy E π of the pion is then imposed by the on-shell condition, while the temporal component ω 1 of q 1 can be continuously varied as indicated in Eq. ( 5), with the temporal component ω 2 of q 2 fixed by energy conservation to ω 2 = E π − ω 1 .The kinematical range accessible on the lattice can be parameterized by the accessible virtualities of the electromagnetic currents, For the Wick rotation to the Euclidean metric to be valid, the virtualities must fall below the two-pion threshold, i.e., q 2 i < 4m 2 π .
For a pion at rest, ⃗ p = 0, it holds that E π = m π and ⃗ q 2 = −⃗ q 1 .In this case, the photon virtualities simplify to The analyticity constraint becomes a constraint on ω 1 , FIG. 3 depicts the kinematical range for 26 different choices of 2 for each of the finite volumes accessible on the three distinct ensembles given in Tab.I.For all three ensembles, the maximum momentum used in the evaluation of the three-point correlation function gives access to virtualities up to |q 2 1,2 | ≈ 1.7 GeV 2 .Using Eqs. ( 1), ( 5) and ( 22) it is straightforward to show that in the rest frame of the pion Ãµν (τ) vanishes when one or both Lorentz indices are temporal, while the spatial components can be written as Inverting the relation gives where Ã(τ) is a scalar under the spatial rotation group.Combining Eqs. ( 1), ( 5) and (25) shows that the form factor can be extracted from a Laplace transform of this scalar quantity as Since the full TFF can be extracted from the scalar Ã(τ) in the rest frame of the pion, we focus on the evaluation of this scalar function for the remainder of this work.The scalar amplitude Ã(τ) is implicitly parameterized by the momentum ⃗ q 1 .All choices of momenta with a fixed value of |⃗ q 1 | 2 can be used to evaluate the form factor as in Eq. ( 26) at identical kinematics.We therefore choose to average the scalar amplitude over such equivalent momenta to increase statistics.

E. Finite-time extent corrections
All ensembles in this work use periodic (anti-periodic) temporal boundaries for the fundamental bosons (fermions).As a result, the pion experiences a periodic boundary condition in time, and the three-point function C µν receives corrections from states propagating backwards around the finite-time extent of our lattices.Denoting the lattice time extent T , including the leading correction gives In the rest frame of the pion, we can apply Bose symmetry and PT symmetry to define a corrected scalar amplitude Ã(τ) (cf.Eq. ( 25)) to account for this leading-order effect, In all further analysis, this corrected form of the lattice amplitude is used as input.

F. Dependence on t π
An example of the scalar amplitude Ã(τ) with the decomposition into contributions from individual Wick contractions is given in FIG. 4. These results, along with those used in the following analysis are based on the choice of interpolating operator time t π = 2.23 fm for ensemble cB211.072.64,t π = 2.18 fm for ensemble cC211.060.80, and t π = 2.28 fm for ensemble cD211.054.96, with the choices made such that t π in physical units is similar on all three ensembles.Note that smaller values of t π lead to a reduction of the statistical errors, however, this requires a careful analysis of excited state effects.To confirm that our choices do not suffer from significant excited state effects, two larger choices of t π were measured on cB211.072.64 and cC211.060.80 and one larger choice on cD211.054.96.We depict an example of the comparison of Ã(τ; t π ) across different choices of t π in FIG. 5.In all cases, Ã(τ; t π ) does not show any systematic dependence on t π across the choices.

G. Tail Fits
The extraction of the transition form factor using Eq. ( 5) in principle requires access to Ã(τ) for arbitrary τ ∈ (−∞, ∞).To estimate the amplitude over the whole temporal range from our lattice data with finite temporal extent, we use two phenomenological models, the Vector Meson Dominance model (VMD) and Lowest Meson Dominance model (LMD) [17,48,49].The VMD and LMD form factors are given, respectively, by and The phenomenological origin of the models suggests particular parameter choices.In the chiral limit and at low energy, the fermionic triangle diagram gives rise to the Adler-Bell-Jackiw (ABJ) anomaly [50,51], which constrains the form factor in these limits to with F π = 92.3(1)MeV the pion decay constant [52], suggesting α = 1/(4π 2 F π ) = 0.274 GeV −1 .Meanwhile, the mass M V = 775 MeV is the ρ-meson mass in the VMD and LMD models, and β = −F π /3 = −0.0308GeV reproduces the leading OPE prediction [27,48,49].However, we treat α, M V and  2 including finite-time extent corrections on the ensembles cC211.060.80 (left) and cD211.054.96(right).We show the full Ã(τ) as black crosses, the connected contribution as orange triangles, and the V-disconnected contributions (multiplied by 10 for easier comparison) with a light, strange and charm vector-current loop as blue, green and yellow circles, respectively.β as free model parameters when fitting the models to our data, using the models as inspiration for the fit form rather than a prediction of the precise asymptotic behavior.
Using the relation between the TFF and amplitude in Eq. ( 1) and inverting the relation in Eq. ( 5), the VMD and LMD form factors can be converted into fit functions for the Euclideantime amplitude Ã as follows.For the LMD fit function one finds in terms of The VMD fit function is obtained by setting β = 0.
The VMD and LMD models should not be expected to fully capture the large-Q 2 behavior of the form factor, or conversely the small-|τ| behavior of Ã(τ).In particular, the Brodsky-Lepage behavior constrains the single-virtual form factor at large Euclidean (space-like) momentum [53][54][55], giving at leading order in α s .This behavior is reproduced by the VMD model, but the LMD model tends to a constant at large Euclidean momenta in the single-virtual form factor.Meanwhile, the operator product expansion (OPE) at short distances [56,57] restricts the double-virtual form factor where both momenta become large at the same time, giving in the chiral limit This behavior is reproduced by the LMD model (with the choice of β given above), but the VMD model in this case falls off as 1/Q 4 [27].Despite these shortcomings at large Q 2 , both models accurately capture the long-distance behavior.As such, we use the lattice data directly for the short-distance region while replacing only the tails with these fits.Details of this procedure are discussed in the following section.

III. TRANSITION FORM FACTOR
To obtain the transition form factors, we need to evaluate the Laplace transform in Eq. ( 26), reproduced here for convenience: where for a given momentum orbit |⃗ q 1 | 2 the squared fourmomenta q 2 1 and q 2 2 are given by Eq. ( 22).The TFF is only directly available on specific sets of orbits q 2 1,2 determined by the kinematics, and it must therefore be interpolated/extrapolated to other points in the (q 2  1 , q 2 2 ) plane.In the following we describe our strategy for the numerical integration in Eq. ( 36) and for the extension of the TFF lattice data to the whole momentum plane using the modified z-expansion.

A. Numerical integration
The tails of Ã(τ) associated with large |τ| are expected to decrease exponentially quickly, with a mass scale predicted by the VMD/LMD models to be E V = M 2 V + |⃗ q 1 | 2 in one tail and E V + m π in the other (see Eq. ( 32)).However, the exponential factor e ω 1 τ enhances one tail of Ã(τ) significantly whenever |ω 1 | is sufficiently large compared to this energy scale.This makes the numerical integration difficult: Since the signal-to-noise ratio decreases when going to larger |τ|, fluctuations in one of the tails due to the small signal-to-noise ratio are exponentially enhanced.We observe that this effect is small for virtualities near the diagonal q 2 1 = q 2 2 , for which ω 1 = m π /2 ≪ E V , but can be severe in the single-virtual cases corresponding to either q 2 1 = 0 or q 2 2 = 0, for which ω 1 = |⃗ q 1 | ≈ E V , especially for the larger momentum orbits.FIG.6 shows a representative comparison of the integrand of the Laplace transform for both the diagonal and single-virtual cases.For choices of virtualities between these two cases, the exponential enhancement gets more and more pronounced when going from the diagonal to the single-virtual case.
These issues are addressed by the introduction of cutoff times τ R cut and τ L cut , with lattice data replaced by a global fit to Ã(τ) for τ ≥ τ R cut and τ ≤ τ L cut .In our analysis, we perform a fit to Eq. ( 32) across data at all available ⃗ q 1 and for values of τ selected from symmetrical fit ranges on both sides of the peak, i.e., τ ∈ [−τ max , −τ min ] ∪ [τ min , τ max ].Variation of the fit ranges allows us to determine a corresponding systematic error.The left and right cutoff times are selected depending on the sign of ω 1 to include as much data as possible on the exponentially suppressed tail of the integrand, in particular by fixing τ L cut = −T/2 for ω 1 ≥ 0 or τ R cut = T/2 for ω 1 < 0. The choice of the cutoff time τ cut is varied to assess the systematic error associated with this choice.
Finally, we filter the choices of cutoff times and kinematics to demand that for a given momentum orbit and value of ω 1 a minimum percentage of F π 0 γγ (q 2  1 , q 2 2 ) must come from the lattice data.This is done to minimize the introduction of a model dependence in the final result.To do so, we define the data content for a given momentum and value of ω 1 to be .
In this work, we always restrict kinematics to only include TFFs for which ∆ latt.≥ 95%.

B. Sampling in the Momentum Plane
As illustrated in FIG. 3, the transition form factor obtained from the lattice is a continuous function of ω 1 for each spatial momentum orbit |⃗ q 1 | 2 .In order to interpolate/extrapolate to the rest of the (q 2  1 , q 2 2 ) plane, we first sample choices of ω 1 on which to evaluate the TFF as the input points for the fit.In this work, we do so by selecting ω 1 corresponding to fixed choices of the ratio q 2 2 /q 2 1 , i.e., by finding the intersection between the available orbits and several diagonal "cuts" through the (q 2  1 , q 2 2 ) plane.Because the underlying lattice data are the same for all choices of ω 1 within each orbit, data at nearby values of ω 1 are strongly correlated.A relatively sparse sampling is therefore possible without sacrificing useful inputs to the z-expansion fits.Both q 2 2 /q 2 1 = 0 and q 2 2 /q 2 1 = 1 are useful choices of ratios to include, the former since the single-virtual transition form factor plays an important role in the evaluation of a π 0 −pole µ , and the latter since we get the best signal-to-noise ratio in the transition form factor there.
We use the following procedure for the sampling in the momentum plane.We determine five values of ω 1 on the highest momentum orbit |⃗ q 1 | 2 = 32(2π/L) 2 , including those corresponding to the diagonal and single-virtual kinematics, such that the arc length of the curve parameterized by ω 1 between neighbouring samples is constant.This fixes five ratios q 2 2 /q 2 1 which are then used to determine the ω 1 on all other spatial momentum orbits.To better cover the region close to single-virtual kinematics, we additionally include a cut q 2 2 /q 2 1 = 0.1, as depicted in FIG. 7. We use several subsets of these cuts in the further analysis in order to determine systematic uncertainties associated with the sampling: all combinations of three cuts including both q 2 2 /q 2 1 ∈ {1, 0} and one of q 2 2 /q 2 1 ∈ {0.88, 0.78, 0.59, 0.1}, and in addition a subset containing the four cuts q 2 2 /q 2 1 ∈ {1.0, 0.88, 0.1, 0.0}.Including FIG. 6. Integrands Ã(τ)e ω 1 τ with |⃗ q 1 | 2 = 6(2π/L) 2 on the ensemble cB211.072.64,comparing diagonal kinematics q 2 1 = q 2 2 (left) and singlevirtual kinematics q 2 2 = 0 (right).We show the integrand from lattice data as black circles and from a correlated LMD fit to Ã(τ) according to Eq. ( 32) as an orange line.The fit was done globally including all momentum orbits (2π/L) 2 ≤ |⃗ q 1 | 2 ≤ 32(2π/L) 2 with fit range [−10, −8] ∪ [8,10], indicated by the grey bands, yielding χ 2 /d.o.f.= 1.01 with correlations taken into account across all orbits and fit range.The red data point indicates the timeslice where the pseudoscalar operator is inserted, while the greyed-out points are those where the time ordering does not satisfy the one from Eq. ( 4).more than four cuts usually leads to an ill-conditioned covariance matrix yielding nonconverging or unstable fits.
For each cut, we impose a threshold for a minimal data content ∆ latt.≥ 95% (see Sec. III A) in the TFF evaluated at each sampled kinematical point, dropping those from further analysis that do not meet the threshold.

C. z-Expansion
We use the modified z-expansion proposed in [28] to interpolate/extrapolate the data from the sampling as described above to the whole kinematical range.This is a modelindependent way of extending the transition form factor to arbitrary photon momenta which is preconditioned to more easily reproduce the form factor structure.Following Ref. [28], the model-independent fit form is constructed by first defining the conformal variables z 1 and z 2 as [58] where Q 2 k = −q 2 k and t c = 4m 2 π indicates the position of the branch cut due to the two-pion threshold.The parameter t 0 can be freely tuned to optimize the rate of convergence.For a given choice of Q 2 max , to best map the region 0 to conformal variables near the origin, the optimal choice of t 0 is which minimizes the maximum value of |z k | in the range [0, Q 2 max ].We use Q 2 max = 4 GeV 2 in the present study, as this is the region yielding the dominant contribution to the integral for a µ (see Sec. IV A).
In terms of these conformal variables, one can expand which is a convergent expansion within the analytic domain |z k | < 1.The coefficients c nm = c mn are symmetric due to the Bose symmetry.By restricting the sum to m, n ≤ N, one can approximate the form factor with accuracy determined by the choice of maximum order N.In addition, the transition form factor can be multiplied by an arbitrary analytic function 2 ) before expanding the resulting product in powers of z k as above.This preconditioning allows the expansion to more easily reproduce the form factor structure.As shown in Ref. [28], the choice where M V ≈ 775 MeV is the ρ-meson mass, leads to a parameterization of the transition form factor which decreases asymptotically as 1/Q 2 in all directions in the momentum plane even at finite values of N, thus satisfying the momentum dependence of Eqs. ( 34) and (35).As shown in Ref. [59], to enforce the appropriate scaling at the two-pion threshold, the expansion should be further modified to fix the derivatives of the transition form factors with respect to z k at z k = −1 to zero.Including the cutoff N, preconditioning, and threshold constraints gives the modified expansion [28] In the following, we either fit the coefficients c nm independently to the data for each ensemble or perform a combined fit to all three ensembles with correction terms for O(a 2 ) lattice artifacts as where we use the cC211.060.In contrast to Refs.[28] and Ref. [29] we perform fully correlated fits with the modified z-expansion.We use N = 1, 2 and all possible combinations of the N = 1 coefficients with one of the N = 2 coefficients.We find that the resulting fits describe the lattice data well, even along the cuts not included in the fit, as seen in the example in FIG. 8 and 9.
Note that the transition form factors obtained from Eq. ( 42) multiplied by Q 2 asymptote to a constant by construction, with predictions for the single-virtual and diagonal case described in Eqs.(34) and (35), respectively.We exclude analyses in which the single-virtual transition form factor multiplied by Q 2 turn out to be asymptotically negative, since this clearly violates those predictions.
We find that correlated z-expansion fits with cutoffs set to either N = 1 or N = 2 already lead to a good quality of fit when including input data from three or four cuts in the momentum plane as described in Sec.III B. For N = 2, at least two of the coefficients c nm are strongly correlated, indicating a degeneracy in the parameterization.Even when including only one of the second-order coefficients in conjunction with all three first-order coefficients, we find slight correlations between some of them.See FIG. 10 for a corner plot depicting correlations in a representative case when including coefficients {c 00 , c 10 , c 11 , c 22 }.
The corresponding correlation matrix is given by +1.00 −0.30 −0.39 +0.31 −0.30 +1.00 −0.05 −0.14 −0.39 −0.05 +1.00 −0.88 +0.31 −0.14 −0.88 +1.00 clearly showing the strong (anti)correlation between c 11 and c 22 , and lesser (anti)correlations between all other coefficients except for c 10 and c 11 which are almost completely uncorrelated.When including all N = 2 coefficients, it turns out that c 20 and c 21 are strongly correlated with c 11 for almost all choices of parameters.To get reliable estimates of the individual coefficients, it is desirable to remove such neardegeneracies in the parameterization by only using a subset of the coefficients which are not strongly correlated.In the final analysis, we therefore determine the preferred values of all quantities studied based on a z-expansion parameterization with minimal degeneracy, but we use variations between these choices of parameterizations to determine an additional systematic error.

IV. OBSERVABLES
The main observable studied in this work is a π 0 −pole µ , i.e., the pion-pole contribution to the muon (g − 2)/2.We further study the two-photon decay width Γ(π → γγ) and the slope parameter b π , both of which are related to F π 0 γγ .
A. Pion-pole contribution to a µ Following Refs.[17] and [18], we use the threedimensional integral representation for the pion-pole contribution where α is the fine structure constant, and 2 ) can be found in FIG. 9.Only TFF data points with at least 95% data content are included (colored in red for the used cuts (q 2 2 /q 2 1 ) ∈ {1.0, 0.88, 0.0}); data in grey with less than 95% data content are shown for illustration only.For Ã(τ) a global LMD fit with fit range [20,21] and τ cut = 23 in lattice units is used.The reduced χ 2 are χ 2 Ã/d.o.f.= 0.86 and χ 2 z−exp./d.o.f.= 1.00.
2 ) in order to highlight the quality of fit at low values of Q 2 2 .The integration is performed over the magnitudes Q 1 , Q 2 of two of the four-momenta and τ = cos θ describing the angle θ between them, with the magnitude of the third fourmomentum fixed by The weight functions w 1 and w 2 are given in Appendix A. Note that w 1,2 are both dimensionless and w 1,2 (Q 1 , Q 2 , τ) → 0 for Q 1,2 → 0 and for τ → ±1.Further, w 2 is symmetric under the exchange of Q 1 and Q 2 .In Ref. [17], the weight functions for the pion are studied and discussed in detail.One finds that the momentum region Q 1,2 ≤ 0.5 GeV is the most important in Eqs. ( 46) and ( 47) for a π 0 −pole µ , which is where we have the strongest data support from our lattice calculation.For two examples illustrating the concentration of the weight functions at low momenta, see FIG. 11.In particular, note that w 2 (Q 1 , Q 2 , τ) is roughly an order of magnitude smaller than By introducing a momentum cutoff Λ, i.e., by replacing in Eqs. ( 46) and ( 47), one can estimate the importance of various momentum regions in the evaluation of Eqs. ( 46) and (47).Though the cutoffs are only directly imposed on Q 1,2 , they imply a cutoff of Q 3 ≤ 2Λ on the third momentum as well.We find that the integrals saturate rapidly, such that we get around 85% of the total result with Λ = 1 GeV and around 90% with Λ = 1.5 GeV, i.e., from momentum regions with strong data support from our lattice calculation.An example of this exercise is shown in FIG.12.

B. Decay Width and Slope Parameter
To leading order in the fine structure constant α, the transition form factors determine the partial decay width through The neutral-pion decay width has been measured in the PrimEx and PrimEx-II experiments [60,61] with a combined result of Γ(π → γγ) = 7.802( 52) stat ( 105) sys (117) tot eV.At our level of statistical precision, uncertainties from radiative corrections to the decay width are negligible and we therefore quote the leading-order result as our estimate of this quantity.Further, the transition form factors can be used to extract the slope parameter thus providing input for determining the electromagnetic interaction radius of the pion.The averaged experimental result for the slope parameter is b π = 1.84(17)GeV −2 [52].These quantities are easily extracted from the form of the z-expansion fit at (q 2  1 , q 2 2 ) = (0, 0).The value of F π 0 γγ (0, 0) comes directly from the z-expansion fit, while the derivative can be acquired by differentiating the form of the z-expansion in Eq. ( 42) with respect to −Q 2  1 , yielding the slope parameter The term in square brackets can be directly evaluated, and the derivative of the conformal parameter is given by

V. MODEL AVERAGING AND ERROR ESTIMATION
We follow different analysis chains from the initial lattice correlation function data to estimates of the quantities a π 0 −pole µ , Γ(π → γγ) and b π .Each such chain is determined by a choice of parameters, namely the choice of fit range and fit model in the fit to Ã(τ), τ cut when constructing the transition form factors and finally the included cuts in the momentum plane in the fit to the modified z-expansion.This yields O(10 3 − 10 4 ) estimates.For the remainder of this section such a choice of parameters is called an "analysis".A priori, only fits with χ 2 Ã/d.o.f.close to 1 are included in an analysis chain for the single ensemble analyses, with 0.84 ≲ χ 2 Ã/d.o.f.≲ 2.05.For for increasing cutoff Λ for ensemble cC211.060.80, using a N = 2 correlated z-expansion fit to (q 2 2 /q 2 1 ) ∈ {1.0, 0.88, 0.1, 0.0}.Only TFF data points where at least 95% of the transition form factors come from lattice data are included, for Ã(τ) a global LMD fit with fit range [7,8] and τ cut = 26 in lattice units is used.The reduced χ 2 are χ 2 Ã/d.o.f.= 0.85 and χ 2 z exp./d.o.f.= 1.02.Indicated is the saturation at Λ ∈ {0.5, 1.0, 1.5} GeV as well as the value of the full contribution for this particular choice of parameters.the combined z-expansion fits the input Ã(τ) fit ranges and τ cut were chosen to be similar in physical units across all ensembles, using t min ∈ [0.48, 1.27] fm, t max ∈ [0.64, 1.43] fm and τ cut ∈ [1.35, 2.07] fm on cB211.072.64 as a reference and remaining within 10% of these physical-units values for the cC211.060.80 and cD211.054.96ensembles.The fits across all ensembles are of good quality, with χ 2 values in the range 1 ≲ χ 2 Ã/d.o.f.≲ 2.4 with an average of χ2 Ã/d.o.f.≈ 1.5.
We calculate the statistical errors for each analysis using the jackknife resampling procedure, finding that there is virtually no autocorrelation between the input data Ã(τ) amongst the available configurations listed in Tab.I. Jackknife resampling is applied over the entire analysis pipeline for each set of analysis choices, ensuring that statistical errors are fully propagated.
To take a weighted average of estimates under various analysis choices in an analysis chain, we use a modified version of the Akaike Information Criterion (AIC) [62,63].Closely following the method introduced in [64] and [65], the model averaging proceeds as follows.For a target observable y, here a π 0 −pole µ , Γ(π → γγ) or b π , we build a histogram from the different analyses, assigning to each analysis a weight given by the AIC.This criterion is derived from the Kullback-Leibler divergence, which measures the distance of a fit function from the true distribution of the points (for a derivation see Ref. [64]).We use the modified AIC introduced in [65], where the χ 2 , the number of fit parameters n par and the number of data points n data describe the fit of interest.The first two terms in the exponent correspond to the standard AIC, and the last term is needed to weigh fits with different lengths in the fit ranges when fitting to Ã(τ) or to weight z-expansion fits with a differing number of cuts in the momentum plane when sampling the transition form factors.
In general, two fitting steps are involved in the determination of each of a π 0 −pole µ , Γ(π 0 → γγ), and b π on a given ensemble: the global fit to Ã(τ) and the z-expansion fit to the transition form factors.A combined weight can be assigned to each combination of fitting choices across both of these steps by multiplying the respective AIC weights wi given in Eq. (53).for the combined fit using c nm = {c 00 , c 10 , c 11 , c 22 }.The orange curve shows the CDF of ≈15'000 different analyses obtained from their AIC weights, the errors on the orange point show the statistical errors for some analyses with a significant weight in the model averaging.This curve corresponds to statistical errors rescaled by λ = 0 while the black curve corresponds to rescaling by λ = 1.The latter gives the central value by the 50% percentile, and the total error as the half variation between the 16% and 84% percentiles shown in gray.To separate the statistical and systematic part of the error, we further calculate the error of the distribution with statistical errors rescaled by λ = 2 (not shown in the figure), as described in the main text.
In the alternative simultaneous continuum and z-expansion fit procedure described in Sec.III C, the model averaging step is instead performed globally across all three ensembles, and therefore a combined AIC weight is assigned by multiplying the weights of the fits to Ã(τ) across all three ensembles as well as the weight of the final global z-expansion fit with lattice artifacts included.In either case, the result is one unnormalized weight wi per analysis.
Given the weights wi , central values m i , and statistical error σ i from each analysis, the global cumulative distribution function (CDF) is expected to be well-described by a weighted combination of normal distribution CDFs.We take the median of the CDF as the central value across all analyses and the half variation between the 16% and 84% percentiles as the total error.To separate statistical and systematic errors, the same procedure can be performed with errors rescaled as σ i → √ λσ i .The variation between λ = 2 and λ = 1 yields a separate statistical and systematic error [65].For an illustration, see FIG. 13.
For a better understanding on the composition of the systematic error, we follow the error budgeting procedure suggested in [65]: For each of the choices made during the analysis chain, e.g., the choice between the VMD or LMD fit to Ã(τ) or the different fit ranges, we first determine the total error for each possible option, varying all other components of the analysis.We then construct a second CDF as above with m i the average of the 16% and 84% percentiles, σ i the total error and w i the sum of the weights coming from this choice.
Using this CDF, we derive the systematic error as described above for the original CDF; this is our result for the systematic error corresponding to the choice.Note that the estimated systematic errors associated with each of the steps of the analysis by this procedure are correlated, thus they do not sum up quadratically to the full systematic error.

VI. LATTICE RESULTS AND CONTINUUM EXTRAPOLATION
Here, our results using the model averaging described in Sec.V are summarized.Comparing to our earlier publication [31] we have refined the analysis by systematically studying different z-expansion fits and by excluding analyses leading to unphysical TFFs as described in Sec.III C.
For the single ensemble analyses, we perform an AIC averaging over different fit models and fit ranges for Ã(τ), different choices of τ cut ≈ [1.3, 2.0] fm on all three ensembles, different samplings in the momentum plane and different choices of z-expansion fits.The results including error budgeting are summarized in Tab.III, and are shown in FIG.16 as diamonds (cD211.054.96),triangles (cC211.060.80), and circles (cB211.072.64),respectively.We note that the error budgets from fit model and fit range are very small, indicating that by restricting the transition form factors in the z-expansion to have at least a 95% contribution from lattice data effectively removes the dependence on the model used to fit the tails of Ã(τ) almost completely.
For the z-expansion fits considered in the combined fitting, the set c nm = {c 00 , c 10 , c 11 , c 22 } with the lattice-artefact correction coefficients δ nm = {δ 00 , δ 10 } leads to the least correlation amongst the coefficients and gives an AIC-averaged fully correlated χ 2 /d.o.f.= 1.36.The AIC-averaged z-expansion coefficients (cf.Eq. ( 42)) in the continuum limit based on the combined fit across all ensembles are given by c 00 c 10 c 11 c 22 0.2220(48) -0.0596(59) -0.050(18) 0.27 (14) with correlation matrix +1.00 −0.46 −0.07 +0.07 −0.46 +1.00 +0.03 −0.08 −0.07 +0.03 +1.00 −0.83 +0.07 −0.08 −0.83 +1.00 and the corresponding corner plot given in FIG.14.This is our preferred continuum result.We note that the coefficients describing the lattice artefacts are very small and compatible with zero, δ 00 = 0.0005 (41) and δ 10 = −0.0085(49).The resulting continuum TFFs for diagonal and single-virtual kinematics are shown in FIG. 15.The figure further indicates the kinematic region directly supported by lattice data, where the fit is expected to be most reliable.For the single-virtual kinematics we find good agreement with the CELLO [19] bins below Q 2 ≲ 2 GeV 2 , while at larger momenta Q 2 ≳ 2 GeV 2 our analysis results tend to be systematically lower than the experimental ones.For the quantities a π 0 −pole µ , Γ(π → γγ) and b π , in order to account for different choices of z-expansion parameters in a more conservative way than by using the AIC, we use a procedure inspired by the method described in Ref. [66].It consists of taking the absolute difference of the central value of the combined fit using the coefficient set c nm = {c 00 , c 10 , c 11 , c 22 } and the weighted average of the results from the other considered fits each weighted with 1/σ 2 tot as an additional systematic error σ sys, z−exp.added in quadrature.We find , Γ(π → γγ) and b π , we could also extrapolate their per-ensemble values in a 2 directly.As depicted in FIG.16, we find that fits linear in a 2 agree well with the combined fit continuum results, albeit with significantly larger errors at a 2 = 0. We note that at the current level of accuracy the lattice artefacts in a π 0 −pole µ are compatible with zero.Hence, a continuum extrapolation with a constant is also possible leading to a better χ 2 /d.o.f. and a significantly smaller error.However, in  order to be on the conservative side, we do not consider this value in our analysis.
Our result for a π 0 −pole µ is compatible with our earlier analysis [31], where we used a continuum extrapolation in the individual ensemble estimates of each observable instead of one coming from a combined fit to the TFF across ensembles.It is also compatible with the recent lattice result a π 0 −pole µ = 57.8± 1.8 stat ± 0.9 sys × 10 −11 from the BMW collaboration in Ref. [29], and a π 0 −pole µ = 59.7 ± 3.6 × 10 −11 from the Mainz group in Ref. [28].Comparing to the dispersive result a π 0 −pole µ = 63.0 +2.7 −2.1 • 10 −11 from Refs.[4,25,26], we find that our result is compatible at the level of 1.6 σ.
data support from lattice CELLO CLEO BaBar Belle FIG. 15.Transition form factor from the fit for diagonal (left) and single-virtual (right) kinematics.For the diagonal kinematics, the OPE prediction for large Q 2 is indicated by the horizontal dashed line while for the single-virtual kinematics, experimental values from CELLO [19], CLEO [20], BaBar [21,22] and Belle [23] are shown.The region with direct support from lattice data is shaded in grey. .Indicated are the statistical and total errors.For comparison, we also show linear fits in a 2 on the per-ensemble data points.In each plot, the points with caps on the errorbars correspond to the combined fit (cross), the cD211.054.96(diamond), the cC211.060.80 (triangle) and the cB211.072.64 (circle) result.In the Γ(π → γγ) and b π plots, the square symbols show the experimental values and their errors, namely Γ(π → γγ) = 7.802(117) eV [61] and b π = 1.84(17)GeV −2 [52].
Finally, for b π , our result is compatible with the experimental result b π = 1.84(17)GeV −2 from Ref. [52] at the level of 1.2 σ.We also agree at the level of 1.6 σ with b π = 1.78(12)GeV −2 from the extraction based on Padé approximants [67] and find a slight tension of 2.1 σ with the dispersive result b π = 1.73(5)GeV −2 from Refs.[25,26].In all cases our central value is higher than these results by O(10%).Further investigation of this quantity from the lattice may therefore be of interest.

VII. COMBINED ANALYSIS OF LATTICE AND EXPERIMENTAL DATA
With our determination of the TFF at the physical point and in the continuum limit we can attempt a combined analysis of our lattice data together with the available experimental data from CELLO [19], CLEO [20], BaBar [21,22] and Belle [23].Such an exercise is interesting for several reasons.Firstly, we can test whether or not the systematic difference, as evident from FIG. 15, between the experimental data for the single-virtual TFF in the momentum region Q 2 ≳ 2 GeV 2 and our prediction based on the extrapolation of the lattice data using the modified z-expansion is indeed significant.To this end we perform a combined fit of the same z-expansion as used for our global fit in Sec.VI, i.e., using the coefficient set c nm = {c 00 , c 10 , c 11 , c 22 } with the corresponding lattice-artefact corrections δ nm = {δ 00 , δ 10 } to the lattice data, while including all the available experimental data up to Q 2 ≃ 35 GeV 2 .In the χ 2 /d.o.f.we use equal weights for each set of experimental data and each set of per-ensemble lattice data.As before, we perform a model average over all analysis chains and obtain the AIC-averaged z-expansion coefficients c 00 c 10 c 11 c 22 0.2277(31) -0.0621(52) -0.1087(85) 0.776 (39) with the correlation matrix given in Appendix B. The resulting TFFs for diagonal and single-virtual kinematics are shown in FIG. 17.We notice the stability and slightly better determination of the z-expansion coefficients c 00 , c 10 and the significant shifts in c 11 , c 22 with substantially reduced errors.This leads to a much more precise determination of the single-virtual TFF at large values of Q 2 ≳ 1.5 GeV 2 which is now fully compatible with the experimental data.The AIC-averaged value of χ 2 /d.o.f.= 1.35 demonstrates that the single-virtual data from the lattice and from experiment are indeed completely consistent with each other and can be described by the same TFF function.For the double-virtual TFF we note the significant difference for Q 2 ≳ 3.0 GeV 2 which is not surprising given the fact that there is no data support in that momentum range.Additional lattice data is needed there in order to further restrict the TFF.
Secondly, we can test the stability of the lattice-driven determination of a (56) While these quantities are now much better determined with errors reduced by factors of 1.4, 1.6 and 3.8, the values are still fully compatible with the lattice-only determinations in Eq. (55) in Sec.VI, with the largest shift observed in a π 0 −pole µ upwards by about 1.5σ.This is easy to understand from the fact that the experimental data shifts the TFF to slightly larger values for Q 2 ≳ 2.5 GeV 2 outside the lattice support, and the fact that this momentum region may still contribute up to 3% to the total value of a π 0 −pole µ , cf., for example, FIG.12.Nevertheless, this exercise shows that the low-energy quantities a π 0 −pole µ , Γ(π → γγ) and b π are not very sensitive to the large-Q 2 behavior of the TFF and can safely be determined from the lattice alone.
Thirdly, for phenomenological analyses where the large-Q 2 behavior of the TFF is important, the lattice-and datadriven parameterization given here is most useful.However, we would like to add the cautionary remark that our analysis does not contain a full systematic analysis of the dependence on higher-orders in the z-expansion fits.We likely expect that the corresponding systematic error leads to an increased total error for the TFFs, in particular for those close to diagonal virtuality where there is no support from lattice or experimental data.

VIII. CONCLUSIONS AND OUTLOOK
We have presented an ab-initio computation of the pion transition form factor at the physical point in the continuum limit using twisted-mass lattice QCD, covering the kinematic range relevant for the extraction of the pion-pole contribution to the HLbL.We are able to include all disconnected Wick contractions contributing to the amplitudes and hence relevant for the calculation of the form factors.
This allows us to provide a parameterization of the transition form factor through the modified z-expansion including the determination of correlations between all coefficients.Our results in the double-virtual and low-Q 2 , single-virtual regimes are complementary to the currently available singlevirtual experimental data.Our data in the latter regime is of particular interest as new experimental data will become available [68] which can be compared with our prediction from first principles.
In the single-virtual regime we find agreement between our calculated pion form factor and the experimental data up to Q 2 2 ≲ 2 GeV 2 , while the extrapolation to Q 2 2 ≳ 2 GeV 2 , based purely on our lattice data, systematically undershoots the experimental data from CELLO [19], CLEO [20], BaBar [21,22] and Belle [23], cf.FIG. 15.In contrast, our results for the partial decay width Γ(π → γγ) and the slope parameter b π are compatible with the experimental values, cf.FIG.16.We note, however, that we do not observe any significant tension between the lattice data and the experimental one, as discussed in Sec.VII.
The main result of this paper, namely the pion-pole contribution to HLbL, is compatible with recent lattice and data-driven results, with a relative total error at thesub-6% level.The main systematic effect not fully accounted for in this work is the effect of working at a fixed finite spatial volume.On general grounds finite-size effects can be estimated in QCD to be of the order O(exp(−m π L)), i.e., they are exponentially suppressed.With our values of m π L = 3.6-3.9we expect them to be at most a few percent.In practice, at the physical lattice sizes of 5.1-5.5 fm employed in this work, finite-size effects appear to be negligible compared to the current precision for a π 0 −pole µ and Γ(π → γγ), as demonstrated in Ref. [29] for similar physical lattice sizes.Since there are no lattice results available from previous work for the slope parameter and its sensitivity to finite-size effects, we can not exclude it to be affected more.We leave the study of finite-size effects as future work, potentially using ETMC ensembles at physical pion mass and the same range of lattice spacings with larger lattice sizes up to ≈ 7.5 fm.
To achieve further precision in future work, several options are available: One can use momentum ⃗ p 0 for the pseudoscalar creation operator (i.e., the moving frame) to give a better coverage of the single-virtual axis.In addition, a fourth physical point ensemble with approximately the same volume as the three used here, but with an even smaller lattice spac- For the diagonal kinematics, the OPE prediction for large Q 2 is indicated by the horizontal dashed line while for the single-virtual kinematics, experimental values from CELLO [19], CLEO [20], BaBar [21,22] and Belle [23] are shown.The region with direct support from lattice data is shaded in grey.
ing, is currently in production and could then also be included in this calculation.

22 FIG. 10 .
FIG. 10.Corner plot for the coefficients c 00 , c 10 , c 11 and c 22 from the correlated z-expansion fit to the TFFs on cC211.060.80 as shown in FIG. 8 and 9 from a jackknife resampling.The corresponding correlation matrix is given in Eq. (44).

22 FIG. 14 .
FIG. 14. Corner plot of the AIC averaged coefficients of the combined correlated z-expansion fit, utilizing bootstrap resampling on each ensemble, where only the coefficients c 00 , c 10 , c 11 and c 22 were used.Notice the anticorrelation between c 11 and c 22 and the less severe anticorrelation between c 00 and c 10 .All other coefficients are virtually uncorrelated.The corresponding correlation matrix is given in Eq. (B1).

TABLE II .
[43]es of Z V and Z A for the ETMC ensembles which were used in this work, as determined in Ref.[43]using a hadronic method.

TABLE III .
Values of measured observables on each ensemble, with systematic error budgeting indicated below each result.

TABLE IV .
Preferred continuum result from the combined fit using c nm = {c 00 , c 10 , c 11 , c 22 } and δ nm = {δ 00 , δ 10 }. σ sys, z−exp.denotes the more conservative systematic error from different z-expansion parameter choices as described in the text.
17G.17.Transition form factor from the combined fit to lattice and experimental data for diagonal (left) and single-virtual (right) kinematics.