On the precision of a data-driven estimate of the pseudoscalar-pole contribution to hadronic light-by-light scattering in the muon g-2

The evaluation of the numerically dominant pseudoscalar-pole contribution to hadronic light-by-light scattering in the muon g-2 involves the pseudoscalar-photon transition form factor F_{P gamma^* gamma^*}(-Q_1^2, -Q_2^2) with P = pi^0, eta, eta^\prime and, in general, two off-shell photons with spacelike momenta Q_{1,2}^2. We determine which regions of photon momenta give the main contribution for hadronic light-by-light scattering. Furthermore, we analyze how the precision of future measurements of the single- and double-virtual form factor impacts the precision of a data-driven estimate of this contribution to hadronic light-by-light scattering.


Introduction
The anomalous magnetic moment of the muon a µ serves as an important test of the Standard Model (SM) [1]. Since several years, there is an intriguing discrepancy of 3 − 4σ between the experimental value [2] and the theoretical SM prediction [1,3]. While this could be a sign of New Physics, the hadronic contributions from vacuum polarization (HVP) and light-by-light scattering (HLbL) have large uncertainties, which make it difficult to interpret the deviation as a clear sign of physics beyond the SM. The hadronic uncertainties need to be reduced and better controlled, also to fully profit from new planned g − 2 experiments [4]. While the HVP contribution [5] can be improved systematically with measurements of σ(e + e − → hadrons), the often used estimates for HLbL a HLbL µ = (105 ± 26) × 10 −11 , [6] (1) a HLbL µ = (116 ± 40) × 10 −11 , [1,7] are both based on model calculations [8][9][10][11][12], 1 which suffer from uncontrollable uncertainties, see also [15]. In this situation, a dispersive approach to HLbL was proposed recently in Refs. [16,17], which tries to relate the presumably numerically dominant contributions from the pseudoscalar-poles and the pion-loop to, in principle, measurable form factors and cross-sections, γ * γ * → π 0 , η, η and γ * γ * → ππ, with on-shell intermediate pseudoscalar states. 2 The hope is that this data-driven estimate for HLbL will allow a 10% precision for these contributions and that the remaining, hopefully smaller contribua e-mail: nyffeler@kph.uni-mainz.de 1 There are attempts ongoing to calculate the HLbL contribution from first principles in Lattice QCD. A first, still incomplete, result was obtained recently [13]. See also the approach proposed in [14]. 2 There have been objections raised at this meeting [18] about the implementation of the dispersive approach for the pion-pole contribution. There should be no form factor at the external vertex [12].

Pseudoscalar-pole contribution
The dominant contribution to HLbL arises, according to most model calculations, from the one-particle intermediate states of the light pseudoscalars π 0 , η, η shown in the Feynman diagrams in Fig. 1. We will evaluate only the pseudoscalar-pole contribution of these two-loop diagrams. In order to simplify the notation, we mainly discuss the neutral pion-pole contribution in this section. The generalization to η and η is straightforward.
3 The weight functions w 1,2 (Q 1 , Q 2 , τ) In Fig. 2 we have plotted the weight functions w 1 (Q 1 , Q 2 , τ) and w 2 (Q 1 , Q 2 , τ) for the light pseudoscalars π 0 , η, η as function of Q 1 and Q 2 for θ = 90 • (τ = 0). Three-dimensional plots for a selection of other values of θ and one-dimensional plots as function of τ = cos θ for some selected values of Q 1 and Q 2 can be found in Ref. [20]. Note that although the weight functions rise very quickly to the maxima in the plots in Fig. 2, the slopes along the two axis and along the diagonal Q 1 = Q 2 actually vanish for both functions [20]. We stress that these weight functions are completely independent of any models for the form factors.
We can immediately see from the weight functions for the pion that the low-momentum region, Q 1,2 ≤ 0.5 GeV, is the most important in the corresponding integrals (7) and (8) for a HLbL;π 0 µ . In The values of the maxima of the weight functions for all light pseudoscalars and the locations of the maxima in the (Q 1 , Q 2 )plane for a selection of θ-values have been collected in Ref. [20]. For w 1 (Q 1 , Q 2 , τ) and θ ≤ 150 • , a ridge develops along the Q 1 direction for Q 2 ∼ 0.18−0.26 GeV (maximum along the line Q 1 = 2 GeV). For larger values of θ, the function looks more and more symmetric in Q 1 and Q 2 . This ridge leads for a constant form factor to an ultraviolet divergence (α/π) 3 C ln 2 (Λ/m µ ) for some momentum cutoff Λ with C = 3(N c /(12π)) 2 (m µ /F π ) 2 = 0.0248, see Refs. [10,11]. Of course, realistic form factors fall off for large momenta and a HLbL;π 0 (1) µ will be convergent. The weight function w 2 (Q 1 , Q 2 , τ) for the pion is about an order of magnitude smaller than w 1 (Q 1 , Q 2 , τ). There is no ridge, since the function is symmetric under Q 1 ↔ Q 2 . For τ near −1, the peak is around Q 1 = Q 2 ∼ 0.14 GeV and broader than the one shown for θ = 90 • . The location Flavour changing and conserving processes The weight functions w 2 (Q 1 , Q 2 , τ) for η and η have a similar shape as for the pion, but the peaks are broader. of the peak moves to lower values Q 1 = Q 2 ∼ 0.04 GeV when τ grows towards +1.
The dependence of the weight functions on the pseudoscalar mass through the propagators shifts the relevant momentum regions (peaks, ridges) in the weight functions, and thus also in HLbL, to higher momenta for η compared to π 0 and even higher for η , see Fig. 2. It also leads to a suppression in the absolute size of the weight functions due to the larger masses in the propagators. This pattern is also visible in the values for the contributions to a µ . For the bulk of the weight functions (maxima, ridges) we observe the following approximate relations (not necessarily at the same values of the momenta and angles) Of course, the ratio of the weight functions is given by the ratio of the propagators and is maximal at zero momenta and at that point equal to the ratio of the squares of the masses, but at zero momenta the weight functions themselves vanish. The combined effect is well described by the relations in Eq. (9). Furthermore, for both η and η , the weight function w 2 is about a factor 20 smaller than w 1 . The peaks for the weight function w 1 (Q 1 , Q 2 , τ) for η and η are less steep, compared to π 0 , and the ridge is quite broad in the Q 2 -direction, so that the weight function is still sizeable, compared to the maximum, for Q 2 = 2 GeV. Furthermore, the ridge falls off only slowly in the Q 1direction. In particular for the η , the ridge for θ ≤ 75 • is almost as big as the maximum out to values of Q 1 = 2 GeV. For w 2 (Q 1 , Q 2 , τ) the peaks are broader and larger momenta contribute, compared to the pion.
For the η-meson, the peak in the weight function The peak for the weight function w 2 (Q 1 , Q 2 , τ) is around Q 1 = Q 2 ∼ 0.14 GeV for τ near −1 as for the pion. The location of the peak moves down to Q 1 = Q 2 ∼ 0.06 GeV when τ is near +1. For the η , the peak in w 1 (Q 1 , Q 2 , τ) now occurs for even higher momenta, The locations of the peaks of w 2 (Q 1 , Q 2 , τ) in the (Q 1 , Q 2 )plane follow a similar pattern as for the η meson.

Relevant momentum regions in a HLbL;P µ
In order to study the impact of different momentum regions on the pseudoscalar-pole contribution, we need, at least for the integral with the weight function w 1 (Q 1 , Q 2 , τ) in Eq. (7), some knowledge on the form factor F Pγ * γ * (−Q 2 1 , −Q 2 2 ), since the integral diverges for a constant form factor. 4 For illustration we take for the pion two simple models to perform the integrals: Lowest Meson Dominance with an additional vector multiplet, LMD+V model, based on the Minimal Hadronic Approximation to large-N c QCD matched to certain QCD shortdistance constraints from the operator product expansion (OPE), see Refs. [10,23] and references therein, and the well-known Vector Meson Dominance (VMD) model. Of course, in the end, the models have to be replaced as much as possible by experimental data on the double-virtual TFF or one can use a DR for the form factor itself [21].
The main difference of the models is a different behavior of the double-virtual form factor for large and equal momenta. The LMD+V model reproduces by construction the OPE, whereas the VMD TFF falls off too fast: Nevertheless for not too large momenta, Q 1 = Q 2 = Q = 0.5 [0.75] GeV, the form factors F π 0 γ * γ * (−Q 2 , −Q 2 ) in the two models differ by only 3% [10%], see [20]. Furthermore, both models give an equally good description of the single-virtual TFF F π 0 γ * γ * (−Q 2 , 0) [23,24]. The LMD+V model was developed in Ref. [23] in the chiral limit and assuming octet symmetry. This is certainly not a good approximation for the more massive η and η mesons. For the η and η meson we therefore simply take the usual VMD model as already done in Refs. [7,10].
The two models yield the following results for the pole-contribution of the light pseudoscalars to HLbL (we only list here the central values) a HLbL;η µ;VMD = 12.5 × 10 −11 .
The results (12) and (13) for the pion-pole contribution in the two models are in the ballpark of many other estimates, but they also differ by 9.4%, relative to the LMD+V result, due to the different high-energy behavior for the double-virtual TFF in Eqs. (10) and (11) for Q 1,2 ≥ 1 GeV. In fact, the pattern of the contributions to a HLbL;π 0 µ is to a large extent determined by the model-independent weight functions w 1,2 (Q 1 , Q 2 , τ), which are concentrated below about 0.5 GeV, up to that ridge in w 1 along the Q 1 direction. As long as realistic form factor models for the double-virtual case fall off at large momenta and do not differ too much at low momenta, we expect similar results for the pion-pole contribution at the level of 15% which is in fact what is seen in the literature [1,10,15]. Nevertheless, due to the ridge-like structure in the weight function w 1 , the high-energy behavior of the form factors is relevant at the precision of 10% one is aiming for.
For the η and η , the results in Eqs. (14) and (15) are as expected from the discussion of the relative size of the weight functions in Eq. (9). The result for η is about a factor 4 smaller than for the pion with VMD. The result for η is only slighly smaller than for η. Note that the normalization of the TFF from Γ(P → γγ) and the momentum dependence due to different vector meson masses for η and η also play a role for the results in Eqs. (14) and (15).
For a more detailed analysis, we integrate in Eqs. (7) and (8)  and display the results, relative to the totals in Eqs. (12)-(15), in Fig. 3.
For the pion, the largest contribution comes from the lowest bin Q 1,2 ≤ 0.25 GeV since a large part of the peaks in the weight functions (for different angles θ) is contained in that bin. More than half of the contribution comes from the four bins with Q 1,2 ≤ 0.5 GeV. In contrast, for the η and η , it is not the bin Q 1,2 ≤ 0.25 GeV which yields the largest contribution, since the maxima of the weight functions are shifted to higher momenta, around 0.3 − 0.5 GeV. Furthermore, more bins up to Q 2 = 2 GeV now contribute at least 1% to the total. This is different from the pattern seen for π 0 . The plots of the weight functions for η and η in Fig. 2 show that now the region 1.5 − 2.5 GeV also is important for the evaluation of the ηand η -pole contributions. The VMD model is, however, known to have a too fast fall-off at large momenta, compared to the OPE. Therefore the size of the contributions a HLbL;η µ and a HLbL;η µ in Eqs. (14) and (15) might be underestimated by the VMD model, which could also affect the relative importance of the higher momentum region in Fig. 3.
Integrating both Q 1 and Q 2 from zero up to some upper momentum cutoff Λ and integrating over all angles θ, one obtains the results shown in Table 1. This amounts to summing up the individual bins shown in Fig. 3.
As one can see in Table 1, for the pion more than half of the final result stems from the region below Λ = 0.5 GeV (59% for LMD+V, 64% for VMD) and the region below Λ = 1 GeV gives the bulk of the total result (86% for LMD+V, 92% for VMD). The small difference between the form factor models for small momenta Q 1,2 ≤ 0.5 GeV is reflected in the small absolute difference for a HLbL;π 0 µ in the two models for Λ ≤ 0.75 GeV. For instance, for Λ = 0.75 GeV, the difference is only 0.8 × 10 −11 , i.e. 1.6%. The faster fall-off of the VMD model at larger momenta beyond 1 GeV, compared to the LMD+V model, leads to a smaller contribution from that region to the total. Therefore we can see in Fig. 3 that the main contributions in the VMD model, relative to the total, are concentrated at lower momenta, compared to the LMD+V model, in particular below 0.75 GeV.
On the other hand, Table 1 shows that for η and η , the region below Λ = 0.25 GeV only gives a small contribution to the total (12% for η, 8% for η ). Up to Λ = 0.5 GeV, we get about half (one third) for η (η ) and the bulk of the results comes from the region below Λ = 1.5 GeV, 96% for the η and 93% for the η -meson.   Table 1. Pseudoscalar-pole contribution a HLbL;P µ × 10 11 , P = π 0 , η, η for different form factor models obtained with a momentum cutoff Λ. In brackets, relative contribution of the total obtained with Λ = 20 GeV. For the calculation of the pseudoscalar-pole contribution a HLbL;P µ with P = π 0 , η, η in Eqs. (7) and (8) the singlevirtual form factor F Pγ * γ * (−Q 2 , 0) and the double-virtual form factor F Pγ * γ * (−Q 2 1 , −Q 2 2 ), both in the spacelike region, enter. We are interested here in the impact of uncertainties of experimental measurements of these form factors on the precision of a HLbL;P µ . See Ref. [19] for a brief overview of the various experimental processes where information on the transition form factors can be obtained. In the following, we will only quote the most relevant and most precise experimental references. Ref. [20] contains a detailed analysis of the current experimental situation.
For the single-virtual form factor F Pγ * γ * (−Q 2 , 0) the following experimental information is available. The normalization of the form factor can be obtained from the decay width Γ(P → γγ) [25,26]. Another important experimental information is the slope of the form factor at the origin. In the spacelike region, the slope and the form factor itself have been measured in the process e + e − → e + e − γ * γ * → e + e − P [27]. The extraction of the slope requires, however, a model dependent extrapolation from rather large momenta Q ∼ 0.5 GeV to Q 2 = 0. The slope and the TFF can also be obtained in the timelike region from the single Dalitz-decay P → γ * γ → + − γ with = e, µ [28]. Of course, for the pion only the decay with an electron pair is possible. For the pion the phase space is, however, rather small and the decay is not very sensitive to form factor effects. The corresponding determinations of the slope and the TFF are rather unprecise. The situation is much better for η and η . However, one then still needs to perform an analytical continuation to obtain the form factor at spacelike momenta. Recently a DR has been proposed in Ref. [21] to determine the single-and double-virtual form factor for the pion. So far, only the single-virtual form factor has been evaluated in this dispersive framework with high precision at low momenta Q ≤ 1 GeV. For the η and η , a dispersive approach for the single-and double-virtual TFF has been presented in Ref. [22].
For the single-virtual form factor we parametrize the measurement errors in Eqs. (7) and (8) as follows: The momentum dependent errors δ 1,P (Q) in different bins are displayed in Table 2, based on the analysis in Ref. [20]. There are currently no experimental data for the form factor available in the spacelike region in the lowest bin 0 ≤ Q < 0.5 GeV for π 0 and η, except from Γ(P → γγ), the slope and timelike data for the η-TFF. For the lowest bin we therefore assume an error, based on "extrapolating" the current data sets and the data that will be available soon from BESIII [29] and maybe from KLOE-2 [24] for the pion. If one uses the DR from Ref. [21] below 1 GeV [30] and the current error on the normalization from the decay width, one obtains (conservatively) the uncertainties in the lowest two bins given in brackets in Table 2. Table 2. Relative error δ 1,P (Q) on the form factor F Pγ * γ * (−Q 2 , 0) for P = π 0 , η, η in different momentum regions. The errors for δ 1,π 0 (Q) and δ 1,η (Q) below 0.5 GeV are based on assumptions.
In brackets for π 0 the uncertainties with a DR for the TFF.
The second ingredient in Eqs. (7) and (8) is the doublevirtual form factor F Pγ * γ * (−Q 2 1 , −Q 2 2 ). Currently, there are no direct experimental measurements available for this form factor at spacelike momenta. For the pion, from the double Dalitz decay π 0 → γ * γ * → e + e − e + e − , one can obtain the double-virtual form factor |F π 0 γ * γ * (q 2 1 , q 2 2 )| at small invariant momenta in the timelike region, but the results are inconclusive [31]. There is indirect information available on the double-virtual TFF from the loop-induced decay P → + − ( = e, µ) [25,32]. Without a form factor at the P − γ * − γ * -vertex, the loop integral is ultraviolet divergent. The relation between this decay and the pseudoscalar-pole contribution a HLbL;P µ has been stressed in Ref. [11] and problems to explain both processes simultaneously with the same model have been pointed out [33].
In this situation, models have been used to describe the form factor F Pγ * γ * (−Q 2 1 , −Q 2 2 ) in the spacelike region and thus all current evaluations of a HLbL;P µ are model dependent. Of course, it would be preferable to replace these model assumptions as much as possible by experimental data. In fact, it is planned to determine the double-virtual form factor, at least for the pion, at BESIII for momenta 0.5 ≤ Q 1,2 ≤ 1.5 GeV and a first analysis is already in progress [29], based on existing data.
In analogy to the single-virtual form factor, we parametrize potential future measurement errors for the double-virtual form factor in Eqs. (7) and (8) in the following, simplifying way: where the assumed momentum dependent errors δ 2,P (Q 1 , Q 2 ) in different bins are shown in Fig. 4. The error estimate for the form factor in each bin has been obtained based on a Monte Carlo (MC) simulation [29] for the BESIII detector using the LMD+V model in the EKHARA event generator [34] for the signal process e + e − → e + e − γ * γ * → e + e − π 0 and the VMD model for the production of η and η . Since the number of events N i in bin number i is proportional to the cross-section σ i (in that bin) and since for the calculation of the cross-section the form factor enters squared, the statistical error on the form factor measurement is given according to Poisson In the lowest momentum bin Q 1,2 ≤ 0.5 GeV, there are no events in the simulation, because of the acceptance of the detector. When both Q 2 1,2 are small, both photons are almost real and the scattered electrons and positrons escape detection along the beam pipe. As a further assumption, we have therefore taken the average of the uncertainties in the three neighboring bins as estimate for the error in that lowest bin. This "extrapolation" from the neighboring bins seems justified, since information along the two axis is (or will soon be) available and the value at the origin is known quite precisely from the decay width. Note that although the form factor for spacelike momenta is rather smooth, it is far from being a constant and some nontrivial extrapolation is needed. For instance, for the pion we get F π 0 γ * γ * (−(0.5 GeV) 2 , −(0.5 GeV) 2 )/F π 0 γ * γ * (0, 0) ≈ 0.5, for both the LMD+V and the VMD model [20].
The MC simulation [29] corresponds to a data sample of approximately half of the data collected at BESIII so far. The simulation included only signal events. Based on a first preliminary analysis of the BESIII data [29] with strong cuts to reduce the background from Bhabha events with additional photons, it seems possible, at least for the pion, that the number of events and the corresponding precision for F π 0 γ * γ * (−Q 2 1 , −Q 2 2 ) shown in Fig. 4 could be achievable with the current data set plus a few more years of data taking. Of course, once experimental data will be available, e.g. event rates in the different momentum bins, there will still be the task to unfold the data to reconstruct the form factor F Pγ * γ * (−Q 2 1 , −Q 2 2 ) without introducing too much model dependence.
Taking the LMD+V and VMD models for illustration, the assumed momentum dependent errors from Table 2 and Fig. 4 impact the precision for the pseudoscalar-pole contributions to HLbL as follows a HLbL;π 0 µ;LMD+V = 62.9 +8.9 −8.2 × 10 −11 While for the pion the absolute variations are different for the two models, as are the central values, the relative uncertainty for both models is around 14%. This will also be visible in the following more detailed analysis. We therefore expect that using other form factor models, and, eventually, using experimental data for the single-and doublevirtual form factors, will not substantially change the following observations and conclusions.  Table 3. Impact of assumed measurement errors δ 1,P (Q) and δ 2,P (Q 1 , Q 2 ) in the form factors F Pγ * γ * (−Q 2 , 0) and F Pγ * γ * (−Q 2 1 , −Q 2 2 ) on the relative precision of the pseudoscalar-pole contributions (first line). Lines 2 − 5 show the effects of uncertainties in different momentum regions below and above 0.5 GeV. Lines 6 − 9 show the impact of potential improvements of some of the assumed errors.