Two and Three Pseudoscalar Production in $e^{+}e^{-}$ annihilation and their contributions to $(g-2)_{\mu}$

A coherent study of $e^{+}e^{-}$ annihilation into two ($\pi^{+}\pi^{-}$,$K^{+}K^{-}$) and three ($\pi^{+}\pi^{-}\pi^{0},\pi^{+}\pi^{-}\eta$) pseudoscalar meson production is carried out within the framework of resonance chiral theory in energy region $E\lesssim 2 \, \mathrm{GeV}$. The work of [L. Y. Dai, J. Portoles, and O. Shekhovtsova, Phys. Rev. D 88,056001 (2013)] is revisited with the latest experimental data and a joint analysis of two pseudoscalar meson production. Hence, we evaluate the lowest order hadronic vacuum polarization contributions of those two and three pseudoscalar processes to the anomalous magnetic moment of the muon. We also estimate some higher-order additions lead by the same hadronic vacuum polarization. Combined with the other contributions from the standard model, the theoretical prediction differs still by $2.7 \,\sigma$ from the experimental value.


I. INTRODUCTION
It is well known that Quantum Chromodynamics (QCD) is successful describing strong interactions. In the high energy region, the correlation functions could be well determined by perturbative QCD. However, the situation becomes more complicated in the low energy region, as the strong coupling constant increases when the energy decreases. including the resonances as new degrees of freedom [3][4][5][6]. The construction of the lagrangian is guided by Lorentz invariance and by chiral, unitary and discrete symmetries, i.e. C-, Pparity conservation. The lack of a coupling that may guide a perturbative expansion in the calculations of the amplitudes, is compensated by a model of the large-N C setting (being N C the number of colours) [7][8][9]. As in χPT, this approach produces the relevant operators in the lagrangian, in terms of Goldstone bosons, resonances and external fields, but leaves undetermined their coupling constants.
One may use experimental data to obtain information on the couplings. However, there is one theoretical tool that has proven efficient in this task: one can extract information on the coupling constants by matching the perturbative Green functions of QCD currents, using the operator product expansion (OPE) at leading order, with those constructed in the RχT framework [10][11][12][13][14][15][16]. Actually, RχT can also match, by construction, with χPT by integrating out the resonances in the Lagrangian [4,17], allowing to relate their coupling constants too.
Low-energy processes with many hadrons in the final state involve final-state interactions (FSI) are notoriously difficult to deal with in a model independent way. The use of dispersive approaches to deal with them is possible in some instances, namely when good phenomenological data are available (see for instance Refs. [30][31][32][33][34] for some recent work). In the framework of RχT, this is also achievable as we did in Ref. [29], where both vector-meson dominance and the anomalous terms were considered in a coherent analysis of the e + e − → π + π − π 0 , π + π − η channels, in the energy region populated by many hadron resonances up to E 2.3 GeV. Here we will revisit that work and will extend it to two pseudoscalar production in the light of the new data.
Recent interest on e + e − annihilation into two and three pesudoscalars is driven by their contribution to the anomalous magnetic moment of the muon a µ = (g µ − 2)/2, with g µ the muon Landé factor. The theoretical prediction of a µ has become a major tour de force in the last years because, on the experimental side, it has been measured with high precision, a exp µ = 11659208.9(5.4)(3.3) × 10 −10 [35,36], and there seems to be a 3.3 σ discrepancy with the standard model (SM) prediction [36]. This fact paves the possibility of bringing out new physics contributions. Within the standard model [37,38], the most important contribution, the electromagnetic, is accurately calculated up to tenth-order α 5 e , a QED µ = 11658471.8951(80)×10 −10 , with very small uncertainty [39]. The electroweak contribution at the two-loop level is also well determined as a EW µ = 15.36(0.1) × 10 −10 [40][41][42][43]. The hadronic contribution is considered as the major source of uncertainty and has two components: hadronic light-by-light scattering (HLBL) and hadronic vacuum polarization (HVP). The HLBL cannot be directly estimated from experimental input, and a combination of different theoretical models has estimated it as a HLBL µ = 10.5(2.6) × 10 −10 [44][45][46]. A comprehensive amplitude analysis on γγ → ππ, KK is done in Refs. [47][48][49][50]. They are indeed the constraints on HLBL where the photons are real. HVP is the largest hadronic contribution and it is related with the cross section of e + e − → anything throughout causality and unitarity.
The present value for the leading order HVP contribution is a HVP,LO µ = 694.0(4.0) × 10 −10 [51]. And the next-to-leading order and next-to-next-to-leading order HVP corrections are derived by considering also higher order hadronic loops, a HVP,NLO µ = −9.87(0.09) × 10 −10 , a HVP,NNLO µ = 1.24(0.01)×10 −10 [52]. The computation of the HVP contribution relies heavily on the available experimental data and, consequently, its improvement will come from the accurate measurement of the electron-positron cross-section.
Comparing the theoretical predictions within the SM and the experimental measurement, there is still a discrepancy, as commented above. There are lots of experimental data available. However, there are discordances among different collaborations, even those with the highest statistics datasets. The study of three pseudoscalar production was carried out in Ref. [29], but recently new experimental measurements of e + e − → π + π − π 0 , π + π − η have become available. SND [53] has given a more precise measurement of e + e − → π + π − π 0 in the energy range 1.05 − 2.00 GeV. BESIII [54] provided a measurement for e + e − → π + π − π 0 in a wide energy range between 0.7 and 3.0 GeV using the Initial State Radiation (ISR) method.
In this paper, we give a coherent analysis of e + e − annihilation into two pseudoscalars π + π − ,K + K − and three pseudoscalars π + π − π 0 , π + π − η based on the former work [29], combined with all the recent experimental measurements. In Sec. II we will briefly update the theoretical framework and give the amplitudes calculated by RχT. In Sec. III, we fit the amplitudes to the experimental data up to 2.3 GeV. In Sec. IV, the leading order HVP contribution to g −2 is estimated. Higher-order hadronic contributions are considered in Sec. V.
Finally, we collect our conclusions in Sec. VI. An appendix collects detailed expressions for the involved form factors and decay widths. II.

A. RχT and further improvements on the form factors
Massless QCD exhibits a chiral symmetry that rules its effective field theory at low energy.
χPT, valid at E M ρ , provides the interaction between the lightest octet of pseudoscalar mesons, and of these with external currents. At higher energies we need to take into account the hadronic resonance states, and a successful phenomenological approach is provided by RχT, which aspects of interest for our case we briefly collect here. We follow the language and notation of Ref. [5].
The structure of the lagrangian has, essentially, three pieces: The first piece involves interaction terms with Goldstone bosons that cannot be generated by integrating out the vector resonance states. They are characterized by a perturbative expansion in terms of momenta (and masses), as in χPT. L V kin involves the kinetic term of the vector resonance states and L V−GB the interaction between Goldstone bosons and vector resonance fields. For the processes that we study in this work only the vector resonance fields will be needed. All of these lagrangians include also external fields coupled to scalar, pseudoscalar, vector, axial-vector or tensor currents. The lowest even-intrinsic-parity O(p 2 ) of the L GB Lagrangian is given by being F the decay constant of the pion and ... indicates the trace in the SU (3) space.
The leading Wess-Zumino-Witten term describing the anomaly with odd-intrinsic-parity is [78,79]. The explicit expression of interest for our work is given by where v σ is the external vector current and Φ the multiplet of Goldstone bosons. Higher orders of the L GB lagrangian will not be considered, as we assume that their couplings are dominated by resonance contributions [80].
The kinetic term of the vector resonance field is given by Here the resonances are collected as SU (3) V octets and have the corresponding properties under chiral transformations. The Lagrangian that involves the interaction between Goldstone bosons and vector resonances, L V−GB , couples the later octets with a chiral tensor constituted by the pseudoscalar nonet and external fields. Hence these chiral tensors obey a chiral counting O(p n ). This allows us to assign a label n to the different pieces as L V... (n) , where the numerator indicates the resonance fields in the interaction terms. We will consider For instance, in the antisymmetric formulation for the spin-one vector resonances that we use, where F V and G V are coupling constants not determined by the symmetry. The rest of terms in Eq. (5) are collected in Ref. [5] for the even-intrinsic-parity terms and Refs. [11,19,29] for those of odd-intrinsic parity. The coupling constants of the interaction terms of L V−GB could be extracted from the phenomenology involving those states. As commented in the introduction the matching between the leading order in the OPE expansion of specific Green functions of QCD and their expressions within RχT is also a useful tool that has been employed in the bibliography [10][11][12][13][14][15][16]. We will implement this procedure as far as it helps in our task. In particular we will use the relations between couplings specified in Ref. [29].
However, the large energy region of study cannot be described fully with only one multiplet of vector resonances V µν . The lightest one is situated around M ρ , i.e. under 1 GeV.
Two other vector multiplets populate the interval 1 GeV E 2 GeV, that we will call V µν and V µν . Their couplings to the pseudoscalar mesons will be defined with respect to the ones of the lightest multiplet as β ππ ,β ππ ,β KK ,β KK , through their poles, as The ρ − ω mixing, required by the e + e − → π + π − process, is reconsidered. While a constant mixing angle δ 0 is enough to describe mixing in the three pseudoscalar case as discussed in Ref. [29]: an energy dependent mixing angle is discussed in Ref. [81], although in the non-relativistic limit and we need to generalize it to the relativistic case. Hence, the energy dependent mixing angle could be parameterized as where |ρ 0 , |ω denote the physical states. Hence the energy dependence of the mixing angle is driven by the resonance propagators. Here M V is the mass of the nonet of vector resonances in the SU (3) limit. We will take M V = M ρ . For the two body final state processes e + e − → π + π − , K + K − , we always take energy dependent mixing mechanism according to Eq. (9). For the three body cases, we adopt two ways. One is to take the same energy dependent ρ − ω mixing mechanism as that of the two body case. This will be Fit I. The other is to use the constant mixing angle δ 0 . This will be our Fit II. Comparison of both fits will unveil the influence of ρ − ω mixing in the analysis of data.

B. Cross sections for two and three pseudoscalars final states
The amplitude for three-meson production in e + e − collisions is driven by the hadronization of the electromagnetic current, in terms of one vector form factor only: being V i µ = qγ µ (λ i /2)q and P = π, η. The Mandelstam variables are defined as s = (Q−p 3 ) 2 , t = (Q − p 1 ) 2 , with Q = p 1 + p 2 + p 3 . The cross section and amplitudes for the three pseudoscalar cases that we are considering, namely e + e − → π + π − π 0 and e + e − → π + π − η, are quite the same as specified in Ref. [29], except for a small change in the treatment of ρ − ω mixing, as illustrated in Sect. II A. The corresponding expressions for the cross-section and the modified form factors for the three pseudoscalar cases are collected in Appendix A.
These form factors depend on several couplings of the RχT lagrangian that are not determined by the symmetry. However, some of them or, at least, relations between them can be established by matching Green functions calculated in this framework with their expressions at leading order OPE expansion of QCD, has been commented before. By implementing these short-distance relations our form factors satisfy both the chiral constraints in the lowenergy region and the asymptotic constraints at the high energy limit (Q 2 → ∞). Hence the only unknown couplings in these form factors will be F V , 2g 4 + g 5 , d 2 , c 3 and α V [29], to be added to the β ππ,KK and β ππ,KK from Eq. (7) and the mixing angles between the octet and singlet pseudoscalar (θ P ) and vector (θ V ) components, defined also in [29].
Two-pseudoscalar final states in e + e − annihilation are given by the corresponding vector form factor with Q = p 1 + p 2 and P = π, K. The energy in the center of mass frame is given by The cross sections σ ππ ≡ σ(e + e − → π + π − ) and σ KK ≡ σ(e + e − → K + K − ) are given by The form factors F π V (Q 2 ) and F K V (Q 2 ) were thoroughly studied in Ref. [82] (see also [83][84][85] for alternative parameterizations) in the case of tau decays. Hence we need to include now the new ρ − ω mixing mechanism, present in e + e − into hadrons. We also extend the described energy region by adding heavier vector multiplets, as commented above. Their expressions are: The functions in Eqs. (13,14) are given by: and Notice that X 1 = 12 sin 2 θ V and X 2 = 12 in the isospin limit. The angles related with the ρ − ω mixing are The Q 2 dependence of resonance widths are a debated issue. A thorough proposal within the chiral framework was proposed in Ref. [86] for wide resonances. We will use this result for Γ ρ (Q 2 ), while a parameterization in terms of the on-shell widths, driven by the two-body phase-space decay will be employed for Γ ρ ,ρ (Q 2 ). The precise expressions are collected in Ref. [29]. Meanwhile the rest of resonances, that are quite narrow, will be taken constant.
Notice that the two-body vector form factors do not include more unknown couplings to those of three-body form factors.

III. COMBINED FIT TO EXPERIMENTAL DATA
As we have seen RχT provides a controlled setting to extract information from experimental data. Part of, but not all, of the couplings have been constrained by demanding that Green functions, in this framework, match the asymptotic behaviour of QCD, within the leading term of the OPE expansion, in the high energy limit. The remaining coupling constants, the mixing angles and resonance masses and on-shell widths are left to be determined from the experimental data of cross sections and widths involving vector resonances.
The unknown couplings include F V , 2g 4 +g 5 , d 2 , c 3 , α V , the phenomenological parameters, β X and β X with X = π, η, ππ, KK, counting for the corresponding strength of the couplings of the V and V [87]. The mixing angles of the pseudoscalar singlet and octet θ P , that of vector singlet and octet θ V , and the ρ − ω mixing angle, the energy dependent δ and/or constant δ 0 are also left free. The masses and widths of resonances belonging to heavier second and third multiplets are also fitted around the central values listed in PDG [36].
The last thirty years of experimental work have been very fruitful getting results for the cross-sections we are interested in, as collected in Sect. I. In order to get results for our parameters we decide to fit the experimental data of cross-sections obtained by dedicated experiments in the last twenty years, i.e. we exclude data older than 2000, with one exception: BESIII [54] measured the cross section of e + e − → π + π − π 0 with high statistics above 1.05 GeV, while it has a relatively large uncertainty below that energy. Thus we do not fit the data points below 1.05 GeV from this dataset. In addition we also fit the PDG figures [36] for the decay widths of vector resonances whose expressions are collected in Appendix A.
Two fits are performed: Fit I uses a uniform energy dependent ρ − ω mixing according to Eq. (9). In Fit II, the two body final state cases take into account the energy dependent ρ − ω mixing, while the other processes are carried out with a constant ρ − ω mixing angle, see Eq. (8). The comparison between cross-sections data and the fit is shown in Figure 1 for the three-pseudoscalar case and Figure 2 for the two-pseudocalar case. The captions in the figures collect all data used in the plots and in the fits.
The global fit includes decay widths of related resonances and their results are shown in Table I. The reported errors are obtained, in quadrature, from two components: one arises from the error propagation of the uncertainty of parameters by MINUIT [88], and the other comes from the statistics with dozens of solutions which could also fit to the experimental data sets well. In general, both Fit I and Fit II provide overall reasonable approximations to the experimental figures quoted in the PDG [36].
The quoted errors in the fitted parameters are provided by the MINUIT code too. In general, the parameters in Fit I and Fit II are consistent with those of Fit 4 in Ref. [29], within a deviation of about 10%. F V , 2g 4 + g 5 , θ V , δ and/or δ are mainly determined by the experimental data under 1.05 GeV, where it has higher statistics and precision. However, the joint fit including the e + e − → K + K − process constrain θ V strongly. This can be understood from the form factor in Eq. (14), because the cross section around the φ peak increases with the descent of θ V . In contrast, the cross section of e + e − → π + π − π 0 around the φ peak decreases when θ V goes down, which could be deduced from the expressions in Appendix A. As a consequence, θ V is about 1 • larger than that of Ref. [29]. The inclusion of e + e − → K + K − process also constrains α V , the higher order correction to the F V coupling arising from SU (3) symmetry breaking. The cross section of e + e − → K + K − increases with rising α V . To confront the theoretical predictions to the experimental data of the cross section of e + e − → K + K − , α V is fixed to be negative. Notice that α V is small as it is higher order correction.
The energy dependent ρ − ω mixing angle δ is determined by the e + e − → π + π − process. From Eq. (13), the cross section of e + e − → π + π − is mainly determined by δ, since F V G V /F 2 = 1 is constrained by the high energy behaviour and θ V could be determined as above. The two mechanisms of ρ − ω mixing adopted in Fit I and II have almost no effects on the three body final state case. There is only a very little difference reflected around the ρ peak in the e + e − → π + π − π 0 process. In the energy region around their masses, ρ and ω mix with a relative phase that results in a larger mode of |F π V | 2 . Hence the magnitude of F V and 2g 4 + g 5 are smaller in Fit I in comparison to Fit II and the results in Ref. [29].
The parameters related with the resonance multiplets are almost the same in Fit I and Fit II, but some of them are different from those of Ref. [29]. They are mainly determined by the energy region above 1.0 GeV. Both e + e − → π + π − and e + e − → ηπ + π − processes are sensitive to the masses and widths of ρ and ρ in this energy region. The e + e − → π + π − data gives relative smaller masses and larger widths of ρ , compared with those provided by the e + e − → ηπ + π − process. Hence the combined fitted ρ mass is about 30 MeV smaller and the ρ width is about 100 MeV larger than those in Ref. [29]. The mass and width of ρ also changes slightly. Consequently, the relative weights of the e + e − → ηπ + π − process β η and β η have sizable changes compared with those in Ref. [29]. Meanwhile, the strengths of the e + e − → π 0 ππ − process β π and β π are similar. Notice that in the two-body processes e + e − → π + π − and e + e − → K + K − , the parameters β ( ) ππ and β ( ) KK turn out to be very small with magnitudes 0.2, as expected by lowest meson dominance [10,[101][102][103][104].
Since d 2 , c 3 and θ P are mainly correlated with the e + e − → ηπ + π − process, they also have sizable changes, while masses and widths of other resonance multiplets are quite the same.
In summary, and as shown in Table II, the fitted masses and widths of heavier multiplets are closer to the experimental average values in PDG [36], due to a combination of updated experimental measurements and the constraints from π + π − and K + K − processes.  Notice however, that the masses and widths of ρ and ρ obtained here correspond to the specific definition of the energy dependent width propagator shown in Eq. (40) of Ref. [29], which may not be used by the experimentalists. Hence a precise comparison with the experimental determinations is not straightforward.
Finally θ V and α V change sizeably with respect to the results of Ref. [29] due to the inclusion of the process e + e − → K + K − , so that the partial widths sensitive to θ V and α V become worse. Nevertheless, these partial widths turn out to be bearable with the experimental data from PDG [36], considering the incertitude associated with the theoret- In general, Fit I and Fit II are almost indistinguishable. There is slight difference shown around the ρ peak in e + e − → π + π − π 0 process at 0.6 < E < 1 (GeV), due to the different parametrization of ρ − ω mixing adopted. Noted that Fit II is a little better in this region, since there is one more parameter δ 0 and the energy dependent mixing mechanism designed for the ππ scattering may not be suitable to be applied in the three pion case, where the three body re-scattering happens. Fit II seems also a little better at the φ peak in the e + e − → π + π − π 0 process. This is because that, F V and 2g 4 + g 5 in Fit II are allowed to have larger values than in Fit I, which can slightly compensate the φ peak in e + e − → π + π − π 0 . As illustrated above, the θ V and α V constrained by the e + e − → K + K − will suppress the φ peak in e + e − → π + π − π 0 . The high energy behaviour of e + e − → π + π − , as shown in Figure 2, is balanced with the e + e − → ηπ + π − process through the mass and width of ρ . where and the kernel function is defined as, with Notice that the lower limit in the integration in Eq. (18) depends on the starting contribution and its O(α e ) order. Hence s thr = m 2 π 0 when including the π 0 γ contribution and s thr = 4m 2 π when starting in the ππ contribution.
It is interesting to notice the 1/s 2 enhancement factor (leading order) of contributions of low energies in a had µ (III). Thus the kernel gives higher weight, in particular, to the lowest lying resonance ρ(770) that couples strong to π + π − . This fact is the reason why the pion pair production e + e − → π + π − gives, by far, the largest contribution to a had µ . However, we are in the position to determine the contributions to the muon anomalous magnetic moment relevant to the three and two pseudoscalar final states that we discussed above. They are shown as a C µ with different energy regions in Table III.
Here a C µ (C = ππ, KK, πππ, ηππ) denotes for the lowest order hadronic vacuum polarization contribution of e + e − → ππ, KK, πππ, ηππ, respectively. The error bars for a C µ are given by the combination of the uncertainty coming from MINUIT [88] and the statistics from dozens of solutions that also fit well the experimental data sets. a µ × 10 −10 Ref. [106] Ref. [29] Ref. [107] Ref. [  It is noted that, although different parameterizations of the ρ − ω mixing are adopted in Fit I and Fit II, the individual contributions of each channel are almost the same. A look back to the Figure 1 shows that the results of Fit I are slightly different from the ones of Fit II around the ρ peak in the e + e − → π + π − π 0 process (see the first three graphs).
However, the total contributions to a πππ µ | ≤1.8 GeV are almost the same, as the contribution of Fit I is slightly larger than that of Fit II on the left hand side of he ρ peak, but it is in the opposite situation on the right hand side of ρ peak. They tend to cancel between each other.
Since there is little difference between the two fits, we will discuss below with Fit II. The a C µ evaluated here are consistent with those in Refs. [51,106,107], within their uncertainty. In addition, a πππ µ | ≤1.8 GeV is also consistent with that evaluated based on the cross section fitted in Ref. [29]. On the other hand, a slightly larger a ηππ µ | ≤1.8 GeV is obtained compared with that of Refs. [29,51]. One has to note that the e + e − → ηπ + π − process has a threshold about 1.2 GeV, and therefore has a larger dependence on the resonance multiplets.
As explained above the largest contribution of the hadronic vacuum polarization comes from e + e − → π + π − . In our theoretical framework, the cross section of e + e − → π + π − below 1 GeV is almost fixed with a small dependence on δ, while other parameters contribute little.
Hence the e + e − → π + π − cross-section shares little uncertainty from the parameters. For the e + e − → K + K − process, only θ V and α V are sensitive, but θ V and α V are in tension with the φ peak in e + e − → π + π − π 0 . Hence, there may be considerable uncertainty of our results.
Since we have fitted up to E = 2.3 GeV, we also listed the corresponding a C µ | ≤2. 3 GeV in Table III. The total contribution is a HVP,LO µ = (699.47 ± 3.38) × 10 −10 (22) from Fit II, in combination with the left channels fitted in Ref. [51]. Note that the four contributions we consider here provide the largest uncertainty among all the channels. Combined with the other contributions (QED [39], EW [40][41][42][43], NLHVP [52], NNLHVP [52], HLBL [44][45][46]) within the SM, we also give an estimation of the anomalous magnetic moment of muon in SM. It is about 5 × 10 −10 larger in total than that in Ref. [51]. Hence our estimation of the discrepancy ∆a µ between the theoretical prediction in SM and that measured by experiment is 0.6σ smaller than that in Ref. [51]. Our estimation of ∆a µ = 20.5 ± 7.7 is 2.7σ smaller than that of the experimental value.

TIONS TO a µ
We can also consider the contribution of the hadronic vacuum polarization to higher-order corrections to the leading result of the previous Sec. IV. These have already been computed in the past at next-to-leading (NLO) order [108] and nex-to-next-to-leading (NNLO) order [52].
In our case, however, we will only consider the contribution of two and three pseudoscalars to HVP, as we have obtained in Sec. III.
NLO contributions correspond to one O(α e ) higher with one and two HVP insertions.
They are given by respectively, where R h (s) has been defined in Eq. (19). The label notation and the kernels K (2a,2b,2c) can be read from Ref. [108]. Notice that the lower limit in the integral is taken 4m 2 π because we are only including the contribution to the cross-section of two and three pseudoscalars.
A higher O(α e ) order with up to three HVP insertions correspond to the NNLO order.
Their contributions can be computed as s, s , s ) .
Our results are shown in Table IV. Since Fit I and Fit II are almost indistinguishable, we would just derive the higher order HVP corrections with Fit II. We also quote the results of Ref. [52], although we remind that the later include all the cross-section and not only the two-and three-pseudoscalar contributions (with √ s ≤ 2.3 GeV) to HVP that we have computed. Hence, the difference between both results can be considered as an estimate of the HVP contributions, that we have not included, and of the higher-energy contribution of the two-and three-pseudoscalar channels. The errors have been estimated in the same way as the leading order contributions to a C µ . It is found that these four processes (with the quoted energy upper limit) account for roughly 70 percent of the higher-order HVP corrections to a had µ .

VI. CONCLUSIONS
Combined with the latest experimental data available for e + e − annihilation into three pseudoscalar cases e + e − → πππ, ππη, we carried out joint fits including the two pseudoscalar cases e + e − → π + π − , K + K − , within the framework of RχT in the energy region up to E 2 GeV. Taking into account the possible different mixing mechanisms of ρ − ω in the three and two pseudoscalar cases, two fits have been performed. In Fit I, we apply a uniform energy dependent ρ − ω mixing parametrization. In Fit II, the energy dependent ρ − ω mixing parametrization is only used in the two pseudoscalar channel, while a constant mixing angle is used in the three body case. Overall very reasonable fits for both cases are found. There is no relevant difference between Fit I and Fit II except for a small difference around the ρ peak in the π + π − π 0 case. This indicates that the ρ − ω mixing mechanism that plays an important role in the two pion case may not be exactly the one to be applied in the three body case. However, it will not affect much the descriptions in the three body case, as well as their contribution to the HVP. Our results have been obtained within a QCD-based phenomenological theory framework with a joint fit of four different channels that restrict mutually each other.
The main hadronic contributions to the muon anomalous magnetic moment come from the lower energy region E < 1.05 GeV of the hadronic vacuum polarization input, where few parameters are dominant. Hence, reliable predictions can be made within our theoretical framework from our previous analyses of the two-and three-pseudoscalar contributions to the e + e − cross-section. Accordingly we have computed the leading-order HVP contribution to the anomalous magnetic moment of the muon by including the four main channels, studied previously, in our estimate. The central value of these four channels to HVP is about 5×10 −10 larger than that of Ref. [51]. In consequence, the discrepancy between SM prediction and the experimental measurement decreases to 2.7σ. As an aside, we have also computed the NLO and NNLO HVP contributions to the anomalous magnetic moment of the muon as given by the two-and three-pseudoscalar contributions to the cross-section. The cross section of the e + e − → π + (p 1 )π − (p 2 )P (p 3 ) process (P a pseudoscalar meson) is driven by the vector form factor in Eq. (10) through being m P = m π , m η , depending on the final state. In Eq. (A.1) the integration limits are: with λ(a, b, c) the Källén's triangle function.
The vector form factors relevant for the e + e − → π + π − π 0 , π + π − η cross-sections are given by: with P = π, η. We give now the expressions for the form factors. When notation is not fully specified we refer to Appendix A.3 of Ref. [29].

b. Three-body decays
The three pion decays of the vector resonances are given by: Finally Ω V is defined by being ε µ V the polarization of the vector meson. Within resonance chiral theory the corresponding reduced amplitudes, Ω V , are: