Pion-pole contribution to hadronic light-by-light scattering in the anomalous magnetic moment of the muon

The $\pi^0$ pole constitutes the lowest-lying singularity of the hadronic light-by-light (HLbL) tensor, and thus provides the leading contribution in a dispersive approach to HLbL scattering in the anomalous magnetic moment of the muon $(g-2)_\mu$. It is unambiguously defined in terms of the doubly-virtual pion transition form factor, which in principle can be accessed in its entirety by experiment. We demonstrate that, in the absence of a direct measurement, the full space-like doubly-virtual form factor can be reconstructed very accurately based on existing data for $e^+e^-\to 3\pi$, $e^+e^-\to e^+e^-\pi^0$, and the $\pi^0\to\gamma\gamma$ decay width. We derive a representation that incorporates all the low-lying singularities of the form factor, matches correctly onto the asymptotic behavior expected from perturbative QCD, and is suitable for the evaluation of the $(g-2)_\mu$ loop integral. The resulting value, $a_\mu^{\pi^0\text{-pole}}=62.6^{+3.0}_{-2.5}\times 10^{-11}$, for the first time, represents a complete data-driven determination of the pion-pole contribution with fully controlled uncertainty estimates. In particular, we show that already improved singly-virtual measurements alone would allow one to further reduce the uncertainty in $a_\mu^{\pi^0\text{-pole}}$.


INTRODUCTION
The anomalous magnetic moment of the muon (g − 2) µ constitutes a highly sensitive probe of physics beyond the Standard Model (BSM). Experimentally, its value is dominated by the final report of the BNL E821 experiment [1], which departs from the SM by about 3σ and thus adds to the intriguing hints for BSM physics observed in the muon sector. Currently, experiment and theory are known at the same level, so that upcoming improved measurements at Fermilab [2] and, potentially, at J-PARC [3] (see also [4]) demand concurrent advances in the SM prediction. If the present discrepancy were to be confirmed, its significance would crucially depend on the theory uncertainties.
In practice, this program requires improvements in two classes of hadronic corrections, see diagrams (a) and (b) in Fig. 1: hadronic vacuum polarization (HVP) and HLbL scattering [5], with QED at five-loop order [6][7][8] (and analytical cross checks at four loops [9][10][11]), electroweak corrections [12,13], and higher-order iterations of HVP [14,15] and HLbL scattering [16] already known sufficiently accurately. As early as [17] it was realized that analyticity and unitarity allow one to express HVP in terms of the cross section σ(e + e − → hadrons), so that improved measurements of the hadronic input directly benefit the SM prediction for (g − 2) µ . Indeed, the most recent compilations [18][19][20] already quote an uncertainty at the same level as or below HLbL scattering, even though tensions particularly in the 2π channel need to be resolved. A similar approach for HLbL scattering, aiming at reconstructing the four-point function from its singularities, has only recently been developed [21][22][23][24][25][26]. The simplest such singularities take the form of singleparticle poles from the exchange of pseudoscalar mesons, with residues uniquely determined by transition form factors that are accessible in experiment, see diagram (c) in Fig. 1. In the context of hadronic models there has been wide agreement that the pion pole emerges as the leading contribution [27][28][29][30][31][32][33][34][35][36][37][38][39][40][41], both based on large-N c arguments and numerical evidence, but its precise definition has been under debate, with variants such as a constant form factor at one vertex [38] or even an off-shell pion [5] being introduced.
In the dispersive reconstruction of the HLbL tensor such ambiguities do not arise, it is the form as derived in [34,35] that reproduces the correct residue (for a vector-meson-dominance transition form factor, this conclusion has also been confirmed using a dispersion relation for the Pauli form factor [42]). In consequence, most recent model-based evaluations have adopted this dispersive definition of pseudoscalar poles [43][44][45], as have calculations of the pion transition form factor within lattice QCD [46] (complementary to lattice-QCD calculations aiming at the full HLbL contribution [47][48][49][50][51]). From an experimental point of view the pion pole can thus be considered similar to HVP, since a measurement of e + e − → e + e − π 0 cross sections for space-like doublyvirtual kinematics fully determines its contribution to (g − 2) µ . In this paper we show that, in the absence of such measurements, existing data for e + e − → 3π, the space-like singly-virtual process e + e − → e + e − π 0 , and the π 0 → γγ decay width already severely constrain the doubly-virtual form factor again by virtue of analyticity and unitarity. This calculation completes a dedicated effort to obtain a fully data-driven determination of the pion-pole contribution to HLbL scattering [52][53][54][55][56][57], including a comprehensive error analysis and matching to short-distance constraints. The conceptual advances presented here will be invaluable to complete a similar program for the η and η ′ poles [58][59][60][61], and thereby all light pseudoscalar single-particle intermediate states in the HLbL contribution to (g − 2) µ .

PION-POLE CONTRIBUTION TO (g − 2)µ
A dispersive approach to HLbL scattering is based on a decomposition of the tensor describing the HLbL process γ * (q 1 )γ * (q 2 ) → γ * (−q 3 )γ(q 4 ) into suitable scalar functions Π i . In terms of these functions the contribution to (g − 2) µ may be written as [26] a HLbL where theΠ i denote linear combinations of the full set of Π i , α = e 2 /(4π), the T i are known kernel functions, and the Euclidean virtualities Q 2 i = −q 2 i are parameterized in terms of a single momentum scale Σ according to [62] TheΠ i are then to be reconstructed in terms of their singularities, most prominently the pole originating from π 0 intermediate states. This singularity only appears in two of theΠ i , and q 2 2 ↔ q 2 3 forΠ 2 , with residues determined by the pion transition form factor F π 0 γ * γ * (q 2 1 , q 2 2 ) (the resulting representation is then equivalent to [35]). This form factor, defined by the matrix element of two electromagnetic currents j µ (x) is thus the central object of this study. Its normalization F πγγ is related to the pion decay constant F π = 92.28(10) MeV [63] by a low-energy theorem [64][65][66] Experimentally, this low-energy theorem has been confirmed at the level of 1.4% in the π 0 → γγ decay [67], which is the uncertainty that we will attach to the central value (5) in this work. We note that this accuracy is likely to improve by a factor of 2 at PrimEx-II [68], although the role of chiral corrections remains to be understood [69]. While there is experimental information for space-like singly-virtual kinematics [70][71][72][73], the available data sets cover primarily large virtualities beyond the low-energy region < ∼ 1 GeV, the latter being most relevant for (g − 2) µ , and doubly-virtual input is lacking completely.
In the following, we construct a form factor representation that proves to be remarkably universal: via the first, dispersive, term, it reproduces all low-energy singularities; the second, small, term parameterizes higher intermediate states that are required to both fully saturate a sum rule for F πγγ and to describe high-energy singly-virtual data; and the third term, based on a pion distribution amplitude, makes the representation respect all asymptotic constraints at O(1/Q 2 ). We begin by demonstrating how the crucial low-energy behavior of F π 0 γ * γ * (q 2 1 , q 2 2 ) can be reconstructed from its 2π and 3π singularities, which in turn can be extracted from e + e − → 2π, 3π.
where q π (s) = s/4 − M 2 π , F V π (s) is the electromagnetic form factor of the pion, and f 1 (s, q 2 ) the partial-wave amplitude for γ * (q)π → ππ. In practice, such an unsubtracted dispersion relation does not fully reproduce the low-energy theorem, the sum rule for F vs (0, 0) = F πγγ /2 being violated by about 10% [55], which can be remedied by invoking a subtracted variant instead. Here, we retain the unsubtracted version given its superior asymptotic behavior, introduce a cutoff s iv in the integral, and restore the sum rule by adding an effective pole as explained below.
The integrand in (8) is reconstructed as follows: F V π (s) is described by an Omnès representation [74] fit to the data from [75]. For the ππ P -wave phase shift we use [76,77], a variant thereof that includes ρ ′ , ρ ′′ in an elastic approximation [53], and the parameterization from [78]. The main complexity, however, is hidden in f 1 (s, q 2 ), whose determination requires the solution of Khuri-Treiman equations [79] for γ * π → ππ [55]. These equations then predict the energy dependence of the partial wave, but the normalization a(q 2 ) needs to be determined from experiment. At q 2 = M 2 ω , M 2 φ this normalization can be extracted from the ω, φ → 3π decay widths, at q 2 = 0 its value is determined by another low-energy theorem [80][81][82][83][84], and for s > 9M 2 π it is accessible in e + e − → 3π. Building upon the parameterization derived in [55], we include a conformal polynomial to better account for inelastic states above 1 GeV, which indeed improves the fit significantly in particular above the φ and allows us to describe the SND [85,86] and low-energy BaBar [87] data with a reduced χ 2 ∼ 1, see Fig. 2. For more details of the representation we refer to [57].
In this way, (8) determines, in principle, the full doubly-virtual form factor, but for the analytic continuation into the space-like region it is advantageous to apply yet another dispersion relation in the isoscalar variable with an isoscalar cutoff s is and an integration threshold s thr = M 2 π 0 (due to the ω → π 0 γ decay). This amounts to a double-spectral representation to describe the low-energy properties of the form factor.
For the singly-virtual asymptotics the situation is more complicated, given that now the low-energy part (10)  contributes and (11) becomes less reliable. However, our representation (10) already saturates the Brodsky-Lepage (BL) limit lim Q 2 →∞ Q 2 F π 0 γ * γ * (−Q 2 , 0) = 2F π at the level of 55%, so that only the remainder needs to be generated by higher intermediate states, which can be conveniently achieved by an effective pole in the doublespectral density. Such a pole amounts to an additional term where the residue g eff is adjusted to restore the sum rule for F πγγ and the mass parameter M eff defines the asymptotic value in the singly-virtual direction without affecting the doubly-virtual behavior at O(1/Q 2 ). The resulting parameters g eff and M eff end up around 10% and 1.5-2 GeV, respectively, consistent with the interpretation as an effective pole that subsumes the effects of higher intermediate states.
For the numerical analysis, we vary g eff within the range defined by the PrimEx experiment [67] and fit M eff to the space-like data [70][71][72][73]. In our formalism, the tension of the BaBar data [72] both with the BL limit and the other data sets manifests itself in a strong sensitivity of the fit results on the lower threshold above which data are included as well as in a deteriorating χ 2 once that threshold is increased. For that reason, we define our central value by the fit to all data sets but BaBar, leading to an asymptotic value almost exactly at the BL limit, but choose an uncertainty band +20 −10 % around the central value, which is sufficiently generous to accommodate a 3σ variation around all fit variants, see Fig. 3, ranging from fits including the BaBar data set to imposing the strict BL limit. Since only data above 5 GeV 2 have been included in the fit, this also defines the prediction for the low-energy region, improving the corresponding result from [55] by the matching to the asymptotic constraints. In particular, we find for the slope parameter The upward shift compared to a π = 30.7(6) × 10 −3 [55] precisely reflects the matching to the asymptotic region, see [57] for details.

RESULTS FOR (g − 2)µ
The final representation (6) then defines the pion-pole contribution to HLbL scattering by means of (1) One major component in the error analysis is the uncertainty in F πγγ , currently known at 1.4% from PrimEx [67], but already the preliminary PrimEx-II update of 0.85% [68] reduces the F πγγ -related and total error to 1.0 and +2.7 −2.1 × 10 −11 , respectively. The uncertainty in the dispersive part (10) is estimated by varying the cutoffs between 1.8 and 2.5 GeV, using several representations for the ππ phase shift and F V π (s) as well as different versions of the conformal polynomial in the e + e − → 3π fit. Next, the uncertainty in the BL limit is estimated by a fit to singly-virtual data as represented in Fig. 3. Finally, the asymptotic piece (13) is evaluated for s m = 1.7(3) GeV 2 , which ensures a smooth matching for q 2 1 = q 2 2 = −Q 2 (note that such a relatively low transition point is indeed expected from light-cone sum rules [97][98][99]).
In contrast to our framework, analytical approaches to the pion transition form factor in the context of (g − 2) µ , such as vector-and lowest-meson dominance [35,43,105], FIG. 4: 2d representation of (Q 2 1 +Q 2 2 )F π 0 γ * γ * (−Q 2 1 , −Q 2 2 ), indicating how the respective asymptotic limits are approached. rational approximants [39,44], or resonance chiral theory [40,45], are restricted to the space-like region, resulting in an interpolation between the on-shell point F πγγ and the space-like singly-virtual data. The main advantage of the dispersive approach concerns the fact that all low-energy singularities are reproduced correctly, which enables us to analyze constraints from all low-energy data in a common framework, in particular time-like data from e + e − → 3π. In this way, the dependence on space-like data for the transition form factor itself is reduced appreciably. Indeed, if one were to take the fit to the complete data base at face value, the BL error would shrink to 0.2 × 10 −11 (with a central value increasing to 63.1×10 −11 ), while the generous band chosen in Fig. 3 allows us to stay completely agnostic of any of the involved systematics without affecting the overall error estimate too severely. In addition, within the dispersive approach we can actually predict the doubly-virtual form factor by virtue of its isospin structure (7) and the double-spectral representation (10), thus removing the systematic uncertainty of having to extrapolate singly-virtual data to doubly-virtual kinematics. Finally, the matching scheme that we have developed in this paper ensures that the asymptotic limits in all directions are correct, see Fig. 4, which is hard to achieve within a resonance model. While our central value lies within the same ballpark as previous calculations, we have therefore, for the first time, provided a robust, comprehensive, and fully data-driven uncertainty estimate.
In general, our error analysis thus shows that, for the pion-pole contribution, the dominant uncertainties are all related to singly-virtual measurements, which points towards opportunities for future improvements. That is, the normalization error in F πγγ is likely to shrink by a factor of 2 at PrimEx-II; the uncertainties in the lowenergy region are currently dominated by the systematics of the representation used for the analytic continuation from the time-like e + e − → 3π data, which could be reduced by including low-energy space-like data, as forthcoming at BESIII [106]; the uncertainties in the BL limit would be reduced substantially if the Belle data for large virtualities were corroborated; and the treatment of the asymptotic region could be improved by comparing to doubly-virtual results from lattice QCD [46]. Irrespective of such future improvements, our result (16) shows that even with currently available data the pion-pole contribution to HLbL scattering in (g − 2) µ is already under very good control, with an uncertainty safely below the level required for the upcoming Fermilab measurement.