Contribution of exclusive $(\pi^0\pi^0, \pi^0\eta, \eta\eta)\gamma$ channels to the leading order HVP of the muon $g-2$

We evaluate the contributions of $(\pi^0\pi^0, \pi^0\eta, \eta\eta)\gamma$ exclusive channels to the dispersion integral of the leading order HVP of the muon anomalous magnetic moment. These channels are included in some way in previous evaluations of the $\pi^0\omega, \eta\omega$ and $\eta\phi$ contributions to $a_{\mu}^{\rm had, LO}$, where the vector resonances (decaying into $\pi^0/\eta+ \gamma$) are assumed to be on-shell. Since the separation of resonance and background contributions in a given observable is, in general, a model-dependent procedure, here we use pseudoscalar mesons and the photon as the $in$ and $out$ states of the $e^+e^- \to (\pi^0\pi^0, \pi^0\eta, \eta\eta)\gamma$ $S$-matrix, such that the cross section contains the interferences among different contributing to the amplitudes. We find $a^{\rm had, LO}_{\mu}(P^0_1P^0_2\gamma)=(1.13\pm 0.13 ) \times 10^{-10}$, where uncertainties stem mainly from vector meson dominance model parameters. Improved experimental studies of these exclusive channels in the whole range below 2 GeV would reduce model-dependency.


INTRODUCTION
During the last two decades, the most accurate measurements of the muon anomalous magnetic moment a µ [1,2] have defied an explanation within the standard model (SM) framework. The reference value of a µ in the SM prediction [3] [2], which is in good agreement with previous results from the BNL 821 collaboration [1]. Forthcoming experimental results from next runs at Fermilab as well as J-PARC [23] and PSI [24] will increase the accuracy reducing the current error by up to a factor of three [3].
The recent measurement of a µ at FNAL [2] arrived simultaneously with a new determination of the hadronic contributions based on lattice QCD [25]. This calculation claims to have reached an accuracy similar to the one of the reference value in the SM (dispersive calculation of the HVP contributions), but it is closer to the experimental value ∆a µ = a exp µ − a SM,LQCD µ = 10.7(6.9) × 10 −10 . Lattice calculations are performed using the fundamental degrees of freedom of QCD to evaluate the HVP contributions; the dispersive evaluations are built up from the sum of cross sections over exclusive hadronic channels to saturate the HVP in the non-perturbative low energy regime. While dispersive calculations of the HVP contributions using the same input data seems to largely agree among them [3,[5][6][7][8][9][10][11], new independent and more precise lattice evaluations may confirm or discard the results of Ref. [25].
If more precise evaluations confirm the difference between lattice and dispersive a had µ results, currently at the 2.1σ level, this will become another interesting anomaly to focus attention on theoretical predictions of a µ . One possible explanation to close the gap may be that some missing or poorly measured low-mass hadronic channels in electron-positron collisions contribute to increase the value of the dispersive integral of the HVP. In this paper we study the contributions of the P 0 1 P 0 2 γ processes (P 1,2 = π or η mesons) to the leading HVP contributions of the muon g − 2 in the SM. These contributions are domi-nated by a rich structure of resonances with masses below 2 GeV. Actually, some of these resonance contributions like ωπ 0 , φη and ωη (with the subsequent radiative decay of vector mesons) have been included in dispersive evaluations of the leading order HVP [9,10], using measurements of the e + e − → V 0 P 0 cross sections (V (P ) will refer hereafter to vector and pseudoscalar mesons). Other exclusive channels involving ω/φ resonances as final states have been reported also in dispersive evaluations of the a LO,had µ [9,10].
Strictly speaking, according to the properties of the S-matrix, the amplitudes involving resonances as incoming/outgoing states are not physical observables [26,27]: only asymptotic physical states n (not resonances) must be included as intermediate states when saturating the unitarity relation: that stems from the S-matrix operator, with S = 1 − iT and SS † = 1 1. This unitarity relation is at the base of the dispersive representation of a had,LO µ and the hadronic cross sections of e + e − annihilations [28][29][30]. Therefore, from a theoretical point of view it is not fully consistent to use resonances as physical final states in hadronic e + e − cross sections, even though it can be a good approximation, particularly for very narrow resonances (see appendix A). This is the main motivation behind the present analysis on P 0 1 P 0 2 γ exclusive channels contributions to a had,LO 1 .
The production cross section of P 0 1 P 0 2 γ states are of the same order in the fine structure constant α as P 0 γ states, with the later being included in evaluations of the HVP contribution (a had,LO µ (π 0 γ + ηγ) 5 × 10 −10 [9,10]). Note that the corresponding non-radiative e + e − → P 0 1 P 0 2 channels are not allowed final states, at least at leading order; therefore, P 0 1 P 0 2 γ do not correspond to their photon inclusive processes. One may think that, given the low threshold for the π 0 π 0 γ its contributions below the 1 GeV region may be enhanced due to the low energy behaviour of the QED kernel in the dispersion integral for a had,LO µ ; however, as it will be shown, the cross sections for P 0 1 P 0 2 γ production is peaked above 1 GeV, leading to suppressed contributions. This property follows from the particular Lorentz structure entering the γ * → P 0 1 P 0 2 γ vertex which leads to e + e − cross sections peaked at center of mass energies above 1.4 GeV. Thus, when those cross sections are inserted into the dispersion integral to evaluate a had,LO µ , the kernel suppression above 1 GeV can be partially 1 Given their large lifetimes compared to hadrons that undergo dominant strong decays, π 0 /η mesons can be considered asymptotic states. Previous calculations of e + e − → π 0 π 0 γ, π 0 ηγ cross sections in the region close to the φ(1020) meson have been provided in Refs. [31][32][33][34]. The corresponding cross sections measurements were reported in [35][36][37], which focus mainly in the hadron mass distribution in φ → P 1 P 2 γ decays. Measurements of the e + e − → π 0 ηγ cross section in the √ s = 1.05 − 2.0 GeV region have been reported by the SND collaboration [38]. More recently, the first measuments of the ηηγ production cross section were reported in [39].
In the absence of experimental data (except for the π 0 ω(→ π 0 γ) channel [40,41]) in the full range below 2.0 GeV, we base our estimate on a Vector Meson Dominance (VMD) model. This model captures the main features of the dynamics of such processes at energies around the resonance regions, and it can be validated with available data as is the case with the measured cross section for e + e − → π 0 π 0 γ [41]. A more sophisticated treatment of the γ * → P 0 1 P 0 2 γ vertex can be done in the framework of resonance chiral theory by including the one-(V P γ) and two-resonances (V V P γ) contributions; Although this analysis is possible it involves a larger set of free parameters associated to the coupling of excited resonances.
We do not consider this and other approaches in the present work.
We organize our paper as follows: after this introduction we describe in Section II the general amplitude and relevant kinematics for the e + e − → P 0 1 P 0 2 γ collisions and introduce some useful notations. In Section III we derive the form factors for the vector-vector contributions to the hadronic vertex. Section IV considers the vector-scalar contributions in the special case of the γ * → π 0 ηγ vertex. In Section V we use available data on the e + e − → π 0 ω(→ π 0 γ) cross section to fit some of the parameters of the model and describe how remaining parameters can be estimated from other data; we also provide our results for the cross sections of different channels. We use the calculated cross sections to compute the dispersive integral and get results for a had,LO µ (P 0 1 P 0 2 γ) in Section VI. Finally, we give our conclusions in Section VII and include two relevant appendices.

II. AMPLITUDE AND KINEMATICS
In S-matrix theory, the quantum amplitudes describe transitions between incoming and outgoing stable states [26]. These initial and final states contain particles which must be described by asymptotic states, i.e. free states that can be defined at times long enough before and after the interaction point. According to this tenet of quantum scattering theory, resonances are not asymptotic states; instead, they appear in transition amplitudes as particles described by propagators of unstable states or poles of the transition amplitudes.
Physical states also form a complete set {|n } of states which satisfy the unitarity condition n |n n| = 1. The unitarity of the S-matrix operator (S = 1 − iT , where T is the transition operator) implies Eq. (1).
Similarly, the use of unitarity in the form of the optical theorem, which allows to relate the HVP of a µ to the cross section for hadron production in electron-positron annihilation via a dispersion relation [28][29][30], requires that only asymptotic states are included in the final states of e + e − annihilations. Experiments have revealed that multihadron production processes are dominated by intermediate resonances which interfere in the squared amplitude.
Owing to interference effects, we can not isolate the observables associated to the production of a given resonance, although it can be a good approximation if the full transition probability is dominated by the production of that resonance [27] (see Appendix A). One such example is precisely e + e − → π 0 π 0 γ, where the intermediate state ω → π 0 γ dominates the cross section.
In this paper, we study how the cross sections behave when one considers the full e + e − → P 0 1 P 0 2 γ processes including all resonances and their interference and we compare our results with the particular case where a single resonance contribution is assumed to dominate the cross section. Our purpose is to reevaluate the HVP contribution to a µ by avoiding the use of resonances as final states.
At the lowest order, the Feynman diagram for this process is depicted in Figure 1. The production amplitude can be presented in the following factorized form: Feynman diagram for e + e − → P 0 1 P 0 2 γ, where P 0 1,2 = π 0 or η . The bubble represents the effects of strong interactions.
where µ =v(p 2 )γ µ u(p 1 ) is the leptonic current and H µ = H µσ * σ = P 0 1 P 0 2 γ|j em µ |0 is the hadronic effective current. The most general form of the hadronic tensor follows from gauge invariance and Lorentz covariance 2 . The form factors A, B, C depend upon the independent Lorentz invariants (q 2 , q 2 , q 2 ) and contains the effects of the strong interactions in the relevant kinematical domain.
The squared amplitude depends upon four independent kinematical variables in addition to q 2 , which is fixed from the total collision energy. Since the P 0 1 P 0 2 γ final states are produced from the s-channel one-photon annihilation of e + e − , the cross section can be written in the following simple form (see for example [42]) where 2 )/(2 q 2 ). The differential cross section in the integrand of Eq. (4) is given by In the following section we consider the VMD model for the hadronic current.
2 This tensor structure is equivalent to the one given in Eq. (2.5) in Ref. [31] with different definition of form factors.
Feynman diagrams describing the γ * → P 0 1 P 0 2 γ vertex in a meson dominance model. Here V (and V ) are intermediate vector meson resonances and S is a scalar meson. Diagrams with exchanged mesons in the final states for diagram (a) must be added to account for Bose statistics (P 0 1 = P 0 2 ) or allowed exchange contributions (P 0 1 = P 0 2 )).

III. FORM FACTORS OF VECTOR-VECTOR CONTRIBUTIONS
In the region In this Lagrangian g is a generic coupling, A µ is the photon field, P (V β ) is the 3×3 matrix of light pseudoscalar (vector) mesons, and Also, one can include isospin and SU(3) breaking effects and try to extract, for lowest lying mesons, the relevant parameters from a global fit to the available data on radiative meson decays, as done for example in Refs. [44] with effective couplings g V P γ (which is related to g defined in (6) for each specifc channel [43,44]). Thus, we can extract the couplings g V P γ from the measured rates of radiative meson decays [43,44] where g V P γ is the coupling for the specific V P γ decay process and P γ the photon threemomentum in the decaying particle's rest frame in each specific decay. The values extracted from the radiative decays of light vector mesons are displayed in the lower part of Table I.
In addition, we need information on the V V P couplings of the radially excited vector mesons V (here V is a light vector meson) and its couplings to photons that enter the γ * → V → V P vertex. Individual measurements of the strong or lepton-pair decays of excited vector mesons needed to determine those couplings, are not reported by the PDG [4].
However, some (model-dependent) analysis of experimental data, mainly from the SND [38,41,45,46], CMD-3 [47], BaBar [48] and BESIII [49] collaborations, allows to extract the ratio of relevant constants g V V P /γ V , where V represents an excited vector meson and em 2 V /γ V its coupling to virtual photons. The relevant product of coupling constants can be extracted from measurements of the cross section at the peak of these resonances which determines the product of their decay rates into V P and lepton pairs [4] through the expression their partial decay widths into the X channel. The values of the X V V P ≡ g V V P /γ V ratios extracted from Eqs. (8) and (20) are given in the upper part of Table I.
With the above ingredients at hand, we can built the amplitude for V V contributions of Figure 2a. The hadronic tensor corresponding to the specific configuration shown in that figure reads: Note that, the last term in Eq. (9) symmetrizes the amplitude for identical mesons in the final state (π 0 π 0 , ηη), and considers the exchange diagram for π 0 η case. In the former case a 1/2! factor must be included in the phase space factor. In the above expression, the form factor F P 1 P 2 γ (u, q 2 , q 2 ; q 2 ) contains information on the production and decay of intermediate resonance states. As expected, the hadronic tensor for vector contributions has the structure derived in Eq. (3).
The squared amplitude for vector contributions will be enhanced at higher c.m.s. energies owing to the quartic momentum dependence Lorentz structure of the hadronic vertex (see Eq. (9)); in addition, it will be further enhanced by the effects of radially excited resonances produced in the s-channel. Owing to this behavior we will include the light and first/second radially excited V resonances in the s-channel, but we keep only the contributions of the lightest vector V resonances decaying into (P 0 2 , P 0 1 )γ final states. Accordingly, we write the form factors for the three processes under consideration as follows (the variables q 2 , q 2 , u and q 2 in the argument of the form factors are omitted): The subindices on the right-hand side refers to the light vector resonances V decaying into (P 0 2 , P 0 1 )γ. An analogous expression to Eq. (11), namely F ηπ 0 γ , must be taken into account for the exchange π 0 ↔ η in the final state. The explicit expressions for each contribution are given in Appendix B.

IV.
FORM FACTORS FOR e + e − → π 0 ηγ In addition to V V contributions discussed in the previous section, scalar resonances can contribute to e + e − → P 0 1 P 0 2 γ as shown in Figure 2b. Among the different final states studied in this paper, only the a 0 (980) and (possibly) its 'excited' scalar states can contribute sizably as intermediate state in the π 0 η channel [35][36][37][38].
The hadronic tensor in this case has a simpler form: This Lorentz structure is in agreement with the general parametrizacion given in Eq. (3).
According to Figure 2b, we need information about the V Sγ and Sηπ interaction couplings. The vertex V Sγ responsible for the scalar resonance production is described by the Feynman rule describing the SP 1 P 2 vertex is given by ig SP 1 P 2 . In terms of these resonances, the form factor for scalar contributions can be parametrized in terms of two scalar resonances a 0 (980), a 0 (1450) (or a 0 , a 0 , respectively, for short) decaying into the π 0 η channel, it becomes (this expression agrees with Eq. (4.1) in Ref. [31] in the case of a single vector and scalar resonance): In the above expressions, One may attempt to extract the relevant couplings of scalar mesons from experimental data. Unfortunately, the experimental information on these decays is rather scarce, if not completely missing 3 . Therefore, we will proceed to use a combination of experimental information, theoretical predictions and make the assumption that only one vector resonance in the s-channel dominates a 0 , a 0 production to provide an estimate of their effects in the cross section: 1. We will assume that the dominant contribution to the e + e − → π 0 ηγ cross section below √ s = 1.2 GeV comes from the γ * → φ(1020) → a 0 (980)[→ π 0 η]γ transition, because both (φ and a 0 ) can be produced on their mass-shell. Therefore, we use the resonance parameters of the a 0 (980) scalar meson as determined in the analysis of the φ → π 0 ηγ hadronic mass distribution measured by KLOE [36] using the resonance model of Ref. 3. We will assume that, in the region of excited vector V resonances, only one of them dominates the γ * → V → a 0 γ vertex. Further, we assume that this vector resonance is the excited state φ = φ(1680). From the following vector-meson dominance relation among the couplings (a similar relation holds for the a 0 γ * γ coupling) and assuming the dominance of the φ (1680), one gets at q 2 = 0, g a 0 γγ = eg a 0 φ γ /γ φ .
Given all the approximations contained in the derivation of scalar couplings, we must take the predicted effects of scalar mesons in the cross section and a had,LO µ (π 0 ηγ) as an indication of their real size.
V. CROSS SECTIONS FOR P 0 1 P 0 2 γ CHANNELS In this section we consider separately the cross sections for the e + e − → (π 0 π 0 , π 0 η, ηη)γ reactions. We focus first, in more detail, on the π 0 π 0 γ channel in order to fix some of the 4 Note that our g a 0 γγ and the one g a 0 γγ used in Ref. [50], are related by g a 0 γγ m 2 a 0 = g a 0 γγ /2. parameters of the model by comparing with available data; this is the channel with the largest cross section among P 0 1 P 0 2 γ final states owing to the large branching ratio for the ω → π 0 γ decay. Thereafter we consider the predictions for the other two channels. given by e + e − → V → π 0 (ρ, ω, φ) → π 0 π 0 γ. Accordingly, the general form of the hadronic tensor was given in Eq. (9) with the specific invariant form factor F π 0 π 0 γ (u, q 2 , q 2 ; q 2 ) = F π 0 π 0 γ ρ + F π 0 π 0 γ ω + F π 0 π 0 γ φ .
The explicit expressions for the different terms are given in Appendix B. Dependence upon the same invariant variables must be understood for each term in the right-hand-side of the above equation.
where m V (Γ V ) denote the mass (width) of resonances. However, following the SND collaboration [40,41] (and our own efforts to achieve a good fit), we use the following expression for the energydependent total width of the ρ(770) meson propagator D ρ (s) = m 2 ρ − s − im ρ Γ ρ (s): where the energy-dependent partial widths are with θ(x) and λ(x, y, z) are the step and Kallen functions, respectively. Although Eq. (17) may look unusual, the opening of new thresholds (like ωπ, KK, · · · ) must be included in the decay width as the invariant mass of the resonance increases.
0.2093 ± 0.0046 [44] TABLE I. Values of model-dependent coupling constants. Entries in the upper part refers to the values extracted from the peak cross sections of e + e − → V → V P as explained in the text using Eqs. (8,20). Values of the middle part are extracted using the VMD expressions for V → P γ and g V P γ couplings (lower part of Table) from Ref. [44].
The form factors given in Appendix B depend upon several parameters: (a) the couplings g V P γ needed to describe (ρ, ω, φ) → P γ decays in the sequence V → P V → P P γ are taken from the fits of Ref. [44], and are listed in the lower part of Table I; (b) the ratio of couplings X V V P ≡ g V V P /γ V shown in the upper part of Table I were extracted from experimental values of the peak cross sections (8) using the theoretical expression (c) the strong V V η couplings for light resonances quoted in the middle part of Table I were extracted by combining the g V ηγ couplings obtained in [44] and the em V /γ V couplings for the γ − V conversion extracted from measured [4] Γ V →e + e − = (4πα 2 /3γ 2 V )m V partial rates and; (d) the masses and decay widths of remaining radially excited vector resonances were taken from [4]. For light vector resonances ρ/ω/φ we assume their masses and widhts world averaged values [4]. In the case of the isovector ρ = ρ(1450) and ρ = ρ(1700) mesons we extract the resonance parameters from a fit to the data of the SND collaboration [41], by assuming their contributions to be complex relative to the lightest ρ(770) vector resonance.
In order to be more explicit, we rewrite the dominant contribution as follows: where X V are taken to be real. The ρ(770) meson propagator D ρ (q 2 ) with the energydependent width is given in Eq. (17).
Parameter SND values [41] This  ) to the e + e − → π 0 ω(→ π 0 γ) cross section data [41], compared to results of Ref. [41] and the PDG values [4]. The † symbol means that the parameter has been kept fixed to their PDG values in the fit.
We can evaluate the cross section using the full S-matrix amplitude by inserting the form factors into Eq. (9), taking into account Bose symmetrization terms, and using Eq. (4) for the cross section. In order to compare to available data from SND [41] on the e + e − → π 0 ω(→ π 0 γ) cross section in the q 2 = 1.05 ∼ 2.0 GeV region 5 , we turn off the first and last terms in Eq. (16). We let as free parameters: the resonance parameters of the ρ(1450), the complex parameters X ρ ωπ , X ρ ωπ (phases φ 1 , φ 2 , respectively) and the g ρωπ coupling.
The third column in Table II collects results of our fit to the cross section data of Ref. [41].
A comparison of the second and third columns in the same Table, shows a good agreement between our results and the ones reported by the SND collaboration [41]. Our fit to the experimental data is shown with a dashed line in Figure 3. In the same Figure,  The amplitude corresponding to V V contributions (Figure 2a) for e + e − → π 0 ηγ must be added with the diagram arising from the exchange of π 0 and η mesons in the final state.
Note however that since the intermediate resonances V and V must be chosen to conserve strong interaction symmetries in the V → ηV and V → π 0 V vertices. Since this exchange contribution does not correspond to the exchange of identical particles in the final states, we do not have to add a 1/2! factor in the phase space.
As it was discussed in Section IV, contributions mediated by scalar mesons can appear in the π 0 η system through the e + e − → γS(→ ηπ) mechanism (S = a 0 , a 0 ). Unfortunately, the situation concerning the experimental information on decay properties of the a 0 (= a 0 (980)) resonance and its nature as a qq, tetraquark or molecular state is not very clear so far [4,[52][53][54][55]. In contrast, the corresponding information for the a 0 (= a 0 (1450)) properties is better known [4].
Despite these limitations, we have attempted an estimate of the effects of scalar resonances. We assume that the dominat contribution is given by the γ * → φ(1020) → a 0 γ FIG. 3. Cross section for the e + e − → π 0 π 0 γ process. The solid line includes all the resonance contributions given in Eq. (16). The dashed line corresponds to the dominance of the ω → π 0 γ decay, the second term in (16). The data points correspond to e + e − → π 0 ω(→ π 0 γ) measured by SND [41].
chain contribution. Similarly, we assume that the a 0 = a 0 (1450) production is dominated by falling of the QED kernel in the dispersive relation.

C. ηηγ final state
The threshold energy for the e + e − → ηηγ is √ s ≈ 1.096 GeV, well above the region of light vector resonances in the s-channel. The form factors for this final state are given in Eq.
(12) and (B1), and the hadronic tensor (9) must include a symmetrization term according to Bose statistics. Using the input couplings shown in Table I, and the convention for the propagators discussed in previous sections, we evaluate the cross section using Eq. (4).
In Figure 5 we plot the e + e − → ηηγ cross section from threshold up to 3.0 GeV. In the absence of experimental information on this decay channel, we assume the different contributions to add coherently with real and positive couplings between different resonance contributions 6 . A dominant peak due to the ρ(1700) is observed and a smaller peak is barely 6 Given that ηηγ contribution to a had,LO µ is at the level of 10 −12 , for now we can keep this approximation FIG. 5. Cross section for the e + e − → ηηγ process. We use the VMD parameters reported in Table I. visible at the φ(2170) resonance position. As expected, the cross section for ηηγ smaller than the one due to π 0 π 0 γ and π 0 ηγ.
The HVP contributions to a µ due to e + e − → P 0 1 P 0 2 γ processes can be written as follows [56,57] a HVP,LO where σ pt is the point cross section for muon-pair production and K(s) is the QED Kernel that can be found, for instance, in Ref. [3].
If we insert the cross sections evaluated in this work into Eq. uncertainty appears in the π 0 π 0 γ contribution. This stems mainly from the uncertainties in the fitted ρ ωπ coupling quoted in Table II [4]. Since we do not use the dataset of all e + e − → π 0 ω(→ π 0 γ) experimental cross sections, our quoted uncertainty for the a had,LO µ (π 0 π 0 γ) channel basically turns out to be larger that the ones quoted by references [9,10] (see discussion below).
We can attempt to make a (risky) comparison with Refs. [9,10], who have provided the evaluations of the π 0 ω(ω → π 0 γ), ηω and ηφ contributions. For the values of a had,LO µ for the latter two channels provided in Refs. [9,10], we add the subsequent decays of (ω, φ) mesons into (π 0 , η)γ decays, which is justified in Appendix A. Under these assumptions we can estimate the P 0 1 P 0 2 γ contributions as follows (B(X) denotes the branching fraction for channel X): Clearly, this represents, at most, an approximation to the complete evaluation. We use the (values obtained in [10] are indicated within square brackets), all numbers in 10 −10 units, and the branching ratios reported in [4] for the radiative decays of vector mesons.
In columns fourth and fifth of Table III we write the values 'estimated' following the above procedure. These values are underestimated with respect to our results and, in the case of the ηηγ channel, by almost one order of magnitude. This is expected since ρη and ρπ 0 exclusive channels are not reported separately and interferences are neglected in Eqs. (23).

VII. CONCLUSIONS
In this paper we have considered the contributions of the neutral e + e − → π 0 π 0 γ, π 0 ηγ and ηηγ exclusive channels to the leading order HVP contributions of the muon anomalous magnetic moment. We evaluate these contributions by considering the full S-matrix amplitude for transitions between these asymptotic states, without cutting intermediate resonances. These decays are not the photon-inclusive channels of e + e − → P 0 1 P 0 2 , P 1,2 = η or π, because such transitions are not allowed (at least at the lowest order in α) and are expected to be of the same order in α as the π 0 γ and ηγ channels. As it is well known [9,10], the later contribute close to 1% to the total contributions of a had,LO µ . We describe the γ * → P 0 1 P 0 2 γ vertex in the framework of Vector Meson Dominance model. We validate this particular model by fitting the available data on the e + e − → π 0 ω(ω → π 0 γ) channel. From the calculated cross sections we evaluate the corresponding dispersion integral and get the following prediction: a had,LO µ (π 0 π 0 γ + π 0 ηγ + ηηγ) = (1.13 +0. 13 −0.14 ) × 10 −10 .
This result is dominated by the π 0 π 0 γ exclusive channel; this is in reasonable good agreement with the evaluation of Refs. [9,10] for the π 0 ω(ω → π 0 γ), where a comparison is possible. The other two contributions are more suppressed and a comparison with existing calculations is not straightforward. Our quoted uncertainty is dominated by errors in the strenght coupling of the ρ → π 0 π 0 γ decay within the VMD model and the particular dataset of e + e − → π 0 ω(→ π 0 γ) measurements [41] used in our analysis.
The cross sections for P 0 1 P 0 2 γ production are peaked in the region populated by excited vector resonances in the s-channel. This introduces important uncertainties in the calculation as long as the information on the parameters and decay properties of excited resonances are rather scarce or not very well known. In order to avoid all the uncertainties related to a particular model, it would be necessary to have better experimental data for these P 0 1 P 0 2 γ final states in electron-positron collisions in the region below 2 GeV.
Of course, this dispersive calculation of a had,LO µ (P 0 1 P 0 2 γ) presented in this paper does not contribute sizably to close the gap with the measured [1,2] and the lattice calculations of Ref. [25]. We address the problem of using exclusive channels with resonances and using them as inputs in the evaluation of a had,LO µ . Using the S-matrix formalism with asymptotic states is important to asses the size of approximations done when one consider resonances as on-shell states and neglects interference with other contributions to the amplitude. It may not be obvious that separating resonance and background contributions from measured observables, is just an approximation. The most clear example that shows that interference effects are important is frequently found in the PDG [4], where the sum over final states involving resonances sometimes exceeds the branching ratios for some specific channels (for example B(D 0 → π + π − π 0 ) = (1.49±0.06)% while i,j B(D 0 → (ρ i (770)π j ) 0 → π + π − π 0 ) = (1.91 ± 0.05)% [4]).

Appendix A: Three Body Scattering Processes
Consider the scattering process A + B → C + D + E. In order to illustrate our point, consider that there are two contributions to the amplitude: M = M R + M B , where the subindex R refers to the production and decays of a resonance R: A + B → C + R(→ D + E) and the subindex B refers to a background (which also may include another less prominent resonance). Accordingly, the cross section contains three terms: σ = σ R + σ B + σ int , where the subdindex 'int' refers to the interference of resonant and background amplitudes.
Isolation of the resonant cross section from the full observable is not possible in general since this requires a good control of background terms. Furthermore, if gauge invariance and gauge-independence is not satisfied by individual contributions in the S-matrix amplitude, the resonance cross section can keep residual gauge-dependence [27]. Note that if background contributions are negigible small in the region around (p C +p D ) 2 ≈ m 2 R ), one can approximate σ(A + B → C + D + E) ≈ σ(A + B → A + R) · BR(R → D + E), where the last factor denotes the branching fraction for the R → D + E decay. This seems to be the case of the π 0 π 0 γ channel discussed in the present paper, where the calculations of the cross section using the full S-matrix for asymptotic states and the aproximation corresponding to σ(e + e − → ω(→ π 0 γ)π 0 ) give very similar results. In this appendix we provide the expressions for the form factors that contribute to the P 0 1 P 0 2 γ transitions as defined in Section III A In the above expressions, the ellipsis in the sum over V s-channel resonances include all possible radial excitations of vector mesons. For identical pseudoscalar mesons in the final state, one needs to exchange q 1 ↔ q 2 in the decay amplitudes, with the corresponding q ↔ q two-particle momenta. Note that for non-identical particles (π 0 η), the form factor for exchanged mesons are not given by the simple exchange of momenta because of the different isospin of π 0 (isovector) and η (isoscalar) mesons.