The η → γ ∗ γ ∗ transition form factor and the hadronic light-by-light η -pole contribution to the muon g − 2 from lattice QCD

,

Introduction.-Radiativetransitions and decays of the neutral pseudoscalar mesons P = π 0 , η and η arise through the axial anomaly and are therefore a crucial probe of the nonperturbative low-energy properties of QCD.The simplest transition to two (virtual) photons, P → γ * γ * , is specified through the transition form factor (TFF) F P →γ * γ * (q 2  1 , q 2 2 ) defined by the matrix element i d 4 x e iq1x 0|T {j µ (x)j ν (0)}|P (q 1 + q 2 ) = µνρσ q ρ 1 q σ 2 F P →γ * γ * (q 2 1 , q 2 2 ), (1) where j µ , j ν are the electromagnetic currents and q 1 , q 2 are the photon momenta.The TFFs determine the partial decay widths to leading order in the fine-structure constant α em through where m P is the pseudoscalar meson mass.Γ(η → γγ) is of particular interest, since it can be used to extract the η − η mixing angles and provides a normalization for many other η partial widths [1].At the same time, there is a long-standing tension between its different experimental determinations through e + e − collisions on the one hand and Primakoff production on the other [2][3][4][5][6][7][8].
The pseudoscalar pole diagrams, depicted in Fig. 1, make the dominant contribution to the HLbL scattering amplitude, with F P →γ * γ * as the key nonperturbative input.Of these diagrams, the π 0 -pole contribution has been estimated based on a dispersive framework [26,51] and on calculations of the pion TFF from lattice QCD [27,52,53] while fits to experimental TFF data using rational approximants have yielded an estimate of all three contributions [24].A preliminary calculation of the η-and η -pole contributions using a coarse lattice was reported in Ref. [54].Experimental results from CELLO [55], CLEO [56], and BaBar [57,58] constrain the spacelike single-virtual 1 GeV 2 , but do not provide data for 0 ≤ Q 2 1GeV 2 or for general double-virtual kinematics.In contrast, these kinematics are the most accessible in lattice QCD.
Evaluating C µν requires the Wick contractions shown in Fig. 3.We evaluate all connected (sub-)diagrams with point-to-all quark propagators.In the P-disconnected diagram, we compute the quark-loop at O † η from propagators based on stochastic volume sources.Unlike in previous lattice QCD studies of the π 0 TFF, here Pdisconnected diagrams of the isospin-singlet η operator are non-zero.The projection onto the η state relies on a delicate cancellation between connected and Pdisconnected diagram contributions, as shown in Fig. 4.
The amplitude Ãµν is then recovered from C µν as where Z η = 0|O η (0, 0)|η( p) is the overlap factor associated with the chosen creation operator.In practice we approximate the limit t η → ∞ by considering three fixed values of t η in the range 0.80 fm t η 1.11 fm.Contamination from excited states and the η meson are suppressed best for the largest value of t η , thus we report the values for Γ(η → γγ), b η , and a η-pole µ from t η ≈ 1.11 fm as the main result and use the remaining choices to check for excited state effects.Statistical noise significantly hinders evaluation of Ãµν (τ ) for large values of |τ |.Following Refs.[27,52,53], we address this issue by performing joint fits of the asymptotic behavior of Ãµν (τ ) for all q 1 to Vector Meson Dominance and Lowest Meson Dominance functional forms [64] with fit windows defined by t i ≤ |τ | ≤ t f .We then integrate over τ as in Eq. ( 5) using numerical integration of the lattice data within the peak region, |τ | ≤ τ c , and analytical integration of the fit form in the tail region, |τ | > τ c .In this work, we consider several choices of τ c in the range 0.16 fm τ c 0.64 fm.
Access to the partial decay width, the slope parameter, and the η-pole HLbL contribution requires an interpolation of the TFF close to the origin and an extrapolation in the quadrant of non-positive photon virtualities.We apply the model-independent expansion in powers of conformal variables advocated in Ref. [27], termed the z-expansion.Analyticity of the form factor below the two-pion thresholds at q 2 1 = 4m 2 π and q 2 2 = 4m 2 π guarantees convergence as the highest power N in the expansion is taken to infinity.Moreover, the expansion is restricted to account for the known threshold scaling and contains preconditioning to more easily capture the expected asymptotic behavior as q 2 1 , q 2 2 → −∞.In practice we find that the N = 2 fit, consisting of six free parameters, already provides a very accurate fit to all lattice results, so we restrict to N ∈ {1, 2} in all subsequent analyses.
The η-pole HLbL contribution has the integral representation [65,66]  with t = cos θ parameterizing the angle between the four-momenta, so that The weight functions w 1 and w 2 are peaked such that contributions to Eq. ( 8) mainly come from the region 0 ≤ Q 1 , Q 2 2 GeV [66].Knowledge of the TFF in the regime of relatively small virtualities is thus sufficient to accurately evaluate a η-pole µ .Finally, we quantify systematic errors associated with the choices of tail-fit model, the parameters (t i , t f ), τ c and the z-expansion order N by the model-averaging procedure detailed in the supplemental material.
Results.-Our lattice results are obtained on the gauge ensemble cB211.072.64 produced by the Extended Twisted Mass Collaboration (ETMC) [67].The seaquark masses for this ensemble are tuned to reproduce the physical charged-pion mass and the strangeand charm-quark masses, with a lattice spacing of a 0.08 fm and a lattice size of L 5.09 fm (m π L 3.62) [63,67].For the valence strange quark we use the mixed action approach in Ref. [61], with valence strangequark mass tuned such that the η meson has physical mass.
We show in Fig. 4 an example of the contributions to C µν (τ, t η ) from the connected and P-disconnected Wick contractions on this ensemble at our largest separation, t η 1.11 fm.The contributions involving strange-quark vector currents are suppressed by a factor ∼ 10 for the connected and ∼ 20 for the P-disconnected contribution compared to those from the light quark vector currents.Contributions from charm-quark vector currents are expected to be even more suppressed [53], as are those from V-disconnected and fully disconnected diagrams (lower two diagrams in Fig. 3).At the presently achievable accuracy these contributions are hence not relevant and are not included in the analysis.
In Fig. 5 we show our results for the TFF as a function of the virtuality in the single-virtual case F η→γ * γ (−Q 2 , 0) (top row) and in the double-virtual case F η→γ * γ * (−Q 2 , −Q 2 ) (bottom row) together with our result from the z-expansion fits.The darker inner band indicates only statistical uncertainties while the lighter outer band includes systematic uncertainties estimated from the variation of fitting choices discussed above.At all virtualities shown, the statistical errors dominate the total uncertainty.In addition to the available experimental data, we also show the Canterbury approximant (CA) result from Ref. [24].We observe reasonable agreement between our results, the experimental data and the CA data.
From the parameterization of the momentum dependence of our TFF data we extract the decay width, slope parameter, and a η-pole µ .As with the estimate of the TFF itself, we repeat the calculation for all choices of the analysis parameters to determine systematic errors associated with tail fits of Ã and the z-expansion.A detailed breakdown is given in the supplemental materials.For the decay width the resulting systematic uncertainty stems mainly from the variation in the fits of the tails of Ãµν (τ ) and τ c , for the slope parameter it is governed solely by the variation of τ c , and for the HLbL pole contribution it is dominated by the conformal fit order N .The total error, however, is always dominated by the statistical uncertainties.We also observe a mild systematic dependence on t η , as detailed in the supplemental materials, which points to the fact that excited-state and possibly η -meson contributions to the transition amplitude are not completely eliminated at the smaller values of t η .We conservatively quote results obtained at our largest value of t η 1.11 fm for which the statistical uncertainty is largest and covers the results at the smaller t η values.
For the leading-order decay width we obtain Γ(η → γγ) = 323(85) stat (22) syst [88] tot eV (9) in comparison to the experimental average 516 (18) eV [1,[3][4][5][6][7].For the slope parameter we find b η = 1.19 (36) stat (16) syst [40] tot GeV −2 (10) to be compared with b η = 1.92(4)GeV −2 from a Padé approximant fit to the experimental results [68] and b η = 1.95(9)GeV −2 from a dispersive calculation [69].Finally, we use the parameterization of our TFF data to perform the integration in Eq. ( 8) and obtain We emphasize that our results are obtained at a fixed lattice spacing and a fixed volume.The present estimates therefore exclude systematic errors associated with finite-volume effects and lattice artifacts.The latter are expected to be of O a 2 Λ 2 QCD with the lattice discretization used here.Lattice artifacts contribute through the bare TFFs, the vector-current renormalization factors (except in b η ) and through the setting of the lattice scale required to convert m µ to lattice units.Both Z V and the lattice scale are determined independently of the quantities considered here [63,67].A quantitative estimate of the lattice artifacts present in a η-pole µ can therefore be obtained by considering the scheme of fixing the renormalization by the physical decay width instead of the hadronic scheme.This gives a η-pole µ;Γ−renorm = 20.9(4.8) stat (2.0) syst • 10 −11 , which differs from a η-pole µ in Eq. ( 11) by 7.7 • 10 −11 , i.e., the effect is similar in size to our total combined statistical and systematic error.
Conclusions and outlook.-Theresults of our lattice QCD calculation of the transition form factor F η→γ * γ * (q 2 1 , q 2 2 ) at physical pion mass have a precision comparable to experimental results in the range where both are available, and demonstrate nice agreement, cf.Fig. 5.For single-virtual kinematics our results provide data at lower photon virtuality than currently accessible by experiments.This includes the region around zero virtuality necessary to study the decay width and slope parameter.Our results in Eqs. ( 9) and ( 10) for these quantities undershoot the experimental (and for b η also theoretical) results by 1.5-2.0standard deviations.
Our lattice computation also provides TFF data for double-virtual (space-like) photon kinematics, which is difficult to access by experiment, making the lattice calculation complementary to the experimental effort.We have made use of this advantage and calculated the η-pole contribution to the anomalous magnetic moment of the muon, a η-pole µ = 13.2(5.2) stat (1.3) syst [5.3] tot • 10 −11 .Our result confirms the currently available data-driven Canterbury approximant estimate [24] and the theoretical model estimates [66,70,71], but does not yet reach the same precision.Nevertheless, it provides important independent support of these estimates.The main shortcoming of our calculation is the use of a single lattice spacing, which will be removed in the future by computations with ETMC gauge ensembles on finer lattices [67,72].
Acknowledgments.We thank Martin Hoferichter, Simon Holz, and Bastian Kubis for helpful discussions.We are also grateful to the authors of Ref. [24] for sharing form factor data produced in their work.This work is supported in part by the Sino- This work [data] This work This work [data] This work [fit] CA FIG. 5. Comparison of the TFF estimated from this work (pink curve) versus the available Fη→γ * γ and Γ(η → γγ) experimental results (blue points) [1,[55][56][57][58] and a Canterbury approximant estimate (cyan curve) [24].For better comparison to features at both small and large Q 2 , the TFFs are plotted both with and without a conventional Q 2 prefactor.

ERROR ESTIMATION AND MODEL AVERAGING
All statistical errors reported in this work are given as 1σ confidence intervals derived from N boot = 2000 bootstrap resamplings of the ensemble of configurations.We find virtually no autocorrelation between the relevant primary data taken on a subset of configurations constituting the ensemble, and the bootstrap bin size is therefore fixed to 1.
During our analysis, we make several choices corresponding to fits of the large-|τ | tails of the amplitude Ãµν (τ ) and of the finite-volume TFF orbits.In particular, the following analysis parameters are varied: 1.The choice between using the Vector Meson Dominance (VMD) or Lowest Meson Dominance (LMD) model to the fit the tail behavior; 2. The window (t i , t f ), determining which regions of the amplitude Ãµν (τ ) are used as inputs to fit the asymptotic tail behavior; 3. The integration cutoff τ c , distinguishing the region |τ | ≤ τ c in which the lattice data is integrated from the region |τ | > τ c in which the analytical tail model is integrated; and 4. The order N of the conformal expansion used to fit the TFFs.
The variation of our estimates with these model choices gives estimates of the systematic errors associated with these steps.We apply the approach of Refs.[74,75] to construct cumulative distribution functions (CDFs) of all final quantities with various subsets of models and with two choices of rescaling parameter λ applied to the systematic error.The various total error estimates, given by the difference between the 16th and 84th percentiles of the CDF in each case, allow an extraction and decomposition of the total uncertainty into statistical, total systematic, and various individual sources.In this approach, weights must be assigned to each model included in the CDF.Weights based on the Akaike Information Criterion [76] derived from χ 2 values of each fit have previously been employed.For the tail of the amplitude, we perform a fit to values of Ãµν (τ ) over sequential choices of τ and across all momentum orbits.For the z-expansion, we perform a fit to values of ) across all orbits at several fixed choices of the ratio Q 2 1 /Q 2 2 .This input data is highly correlated, and determining the correlated χ 2 therefore requires a very precise estimate of nearly degenerate covariance matrices of both the tail fits and z-expansion fits.Even for fits to small windows (t i , t f ) and few choices of orbits, we did not find estimates of the χ 2 values to be accurate in our preliminary investigations.Instead, in this work we perform fits using much more stable uncorrelated fits and make the conservative choice to use a uniform weighting of all possible models in the CDF method.This can be expected to overestimate the systematic error associated with model variation.
The decomposition of uncertainties is detailed in Table I for all three final physical quantities studied in this work.Due to correlations between the total error estimates in each case, the decomposition does not simply add in quadrature, but nevertheless gives an estimate of which components of the error dominate the error budget.Unsurprisingly, the dominant sources of systematic errors vary depending on the observable considered.For the η-pole contribution to the HLbL, the biggest source of systematic error is the conformal fit used to extrapolate the  For reference, the values are respectively compared against estimates from the PDG [1], Padé approximant (PA) fits to experimental data [68], and the VMD model [66] and Canterbury approximant (CA) experimental fits [24].
2 ) from the low-virtuality orbits accessible on the lattice to the full plane of spacelike (q 2 1 , q 2 2 ).This indicates that, despite the important contributions to a η-pole µ from low virtualities, the large uncertainties in the nearly unconstrained higher virtualities can still affect the estimate of a η-pole µ from lattice data alone.Incorporating some information about asymptotic scaling of the TFF at large virtualities is therefore an interesting prospect for future work.The other two quantities, Γ(η → γγ) and b η are directly related to the behavior of the TFF at q 2 1 = q 2 2 = 0 and can therefore be seen to be much less sensitive to the conformal extrapolation.Instead, here the choices used to fit the tails of the amplitude Ãµν (τ ) dominate the systematic errors.Nonetheless, we find that the uncertainties in all three quantities are almost entirely given by the statistical error, which always far outweighs the systematic errors.

DEPENDENCE ON tη
In Fig. S1 we show the dependence of the partial decay width Γ(η → γγ), the slope parameter b η , and the η-pole contribution a η-pole µ on the choice of t η which denotes the imaginary time location of the creation operator O † η (−t η ) for the η meson, to be compared with imaginary time coordinates of the currents j µ (τ ) and j ν (0).The outer error bar denotes the total error, while the inner one shows the statistical error only.It is clear that the total error is dominated by the statistical one in all cases and for all t η considered in this calculation.For all three quantities we observe a mild systematic trend with growing t η which may be an indication that excited state and η -meson contributions to the transition amplitude, and hence to the quantities shown here, may still be present at the smaller values of t η .Since we are interested in the limit t η → ∞ we conservatively quote the results for the largest available t η for which the statistical error is largest and covers the results at the smaller values of t η .

INTERPOLATION OF THE η STATE
The η-meson state is the lowest-lying eigenstate of the twisted mass lattice Hamiltonian in the channel with quantum numbers I G J P C = 0 + (0 −+ ).The exact interpolating field to project onto the η eigenstate in the lattice calculation is unknown.However, it is sufficient that it can be written as a linear combination of the quark-model octet-and singlet-pseudoscalar operators Nevertheless, the η-meson state is the unique ground state of lowest mass, and propagation in Euclidean time systematically suppresses the contribution of the η -meson and excited states lying higher in the spectrum.This suppression scales exponentially as exp (−(M − m η )t), in terms of the Euclidean time evolution t and the relative energy gap between the mass M of the higher state and m η .This applies to all two-and three-point correlation functions used in this work.Thus for sufficiently long Euclidean time propagation, the projection onto the η-meson state is achieved by our choice in Eq. ( 6) of O † 8 as the creation operator for the three-point function.

FIG. 3 .
FIG.3.Wick contractions contributing to Cµν (τ, tη).The second connected diagram with quark propagators running in the opposite direction and the second V-disconnected diagram with a loop at jµ(τ ) are omitted for brevity.
German collaborative research center CRC 110 and the Swiss National Science Foundation (SNSF) through grant No. 200021 175761, 200021 208222, and 200020 200424.We gratefully acknowledge computing time granted on Piz Daint at

TABLE I .
Decomposition of uncertainties in the reported values of the three quantities studied at the single lattice spacing and volume used in this work.The results and uncertainties are based on the conservative choice tη/a = 14 corresponding to tη = 1.11 fm.