Three-pion contribution to hadronic vacuum polarization

We address the contribution of the $3\pi$ channel to hadronic vacuum polarization (HVP) using a dispersive representation of the $e^+e^-\to 3\pi$ amplitude. This channel gives the second-largest individual contribution to the total HVP integral in the anomalous magnetic moment of the muon $(g-2)_\mu$, both to its absolute value and uncertainty. It is largely dominated by the narrow resonances $\omega$ and $\phi$, but not to the extent that the off-peak regions were negligible, so that at the level of accuracy relevant for $(g-2)_\mu$ an analysis of the available data as model independent as possible becomes critical. Here, we provide such an analysis based on a global fit function using analyticity and unitarity of the underlying $\gamma^*\to3\pi$ amplitude and its normalization from a chiral low-energy theorem, which, in particular, allows us to check the internal consistency of the various $e^+e^-\to 3\pi$ data sets. Overall, we obtain $a_\mu^{3\pi}|_{\leq 1.8\,\text{GeV}}=46.2(6)(6)\times 10^{-10}$ as our best estimate for the total $3\pi$ contribution consistent with all (low-energy) constraints from QCD. In combination with a recent dispersive analysis imposing the same constraints on the $2\pi$ channel below $1\,\text{GeV}$, this covers nearly $80\%$ of the total HVP contribution, leading to $a_\mu^\text{HVP}=692.3(3.3)\times 10^{-10}$ when the remainder is taken from the literature, and thus reaffirming the $(g-2)_\mu$ anomaly at the level of at least $3.4\sigma$. As side products, we find for the vacuum-polarization-subtracted masses $M_\omega=782.63(3)(1)\,\text{MeV}$ and $M_\phi=1019.20(2)(1)\,\text{MeV}$, confirming the tension to the $\omega$ mass as extracted from the $2\pi$ channel.


Introduction
Three-particle decays subject to strong final-state interactions are notoriously difficult to describe in a fully model-independent way, i.e., without assumptions on intermediate states of the decay or other approximations of the hadron dynamics. One of the simplest examples is the three-pion decay of vector mesons, V = ω, φ, which phenomenologically is dominated by the ρ(770) resonance formed in the final-state rescattering of the pions. However, a description beyond a simple isobar model is challenging, especially given that the decay is out of reach for low-energy effective field theories. A strategy to control the pion final-state interactions based on analyticity and unitarity was first developed in the context of K → 3π [1] and applied to ω → 3π as early as [2]. These Khuri-Treiman (KT) equations have since become a standard tool in three-particle decays, with recent applications specifically to ω, φ → 3π decays in [3][4][5], in part triggered by significant progress in the determination of the ππ phase shifts that are required as crucial input in the solution [6,7]. A detailed, model-independent understanding of hadronic amplitudes can have significant impact beyond low-energy QCD itself, most notably in low-energy searches for physics beyond the Standard Model (SM) such as the anomalous magnetic moment of the muon a µ = (g − 2) µ /2, whose SM prediction currently disagrees with experiment [8] a exp µ = 116 592 089(63) × 10 −11 (1.1) at the level of 3-4σ. To confront the SM with upcoming experiments at Fermilab [9] and J-PARC [10], one needs to be able to control the theoretical uncertainties at a commensurate level. In this context, the issue of hadronic modeling is most severe in the hadronic light-by-light (HLbL) contribution, for which a data-driven dispersive approach has only been recently developed [11][12][13][14][15][16][17][18]. 1 In contrast, the leading hadronic contribution, hadronic vacuum polarization (HVP), is, in principle, fully determined by the cross section for e + e − → hadrons [22,23], and indeed a combination of the analysis of exclusive channels, inclusive data, and perturbative-QCD constraints are used for current estimates of the HVP contribution [24][25][26][27]. 2 However, only the compilations from [24,25] are exclusively based on the direct integration of the data, while [26,27] do involve some model assumptions, in particular for the ω and φ contributions. In general, tensions among data sets are typically taken into account by a local error inflation. For the lowest-multiplicity channels that dominate HVP at low energies the available constraints from analyticity and unitarity (as well as low-energy theorems) are powerful enough to define a global fit function that the data need to follow if consistent with all QCD constraints. Such an approach to the 2π channel below 1 GeV has recently been completed in [36], relying on a close relation between ππ scattering, the pion vector form factor, and the HVP integral [37][38][39][40][41]. As a result, it was found that despite the known tension between BaBar [42,43] and KLOE [44][45][46][47], each data set by itself is consistent with QCD constraints, and a global fit then defines an average that only uses as additional input information on the covariance matrices as provided by experiment.
Here, we extend this strategy to the 3π channel, which produces both the secondlargest contribution to the total HVP value and its uncertainty. Instead of the pion vector form factor, the underlying hadronic amplitude becomes γ * → 3π, which was studied in detail in the context of the pion-pole contribution to HLbL scattering [17,18,[48][49][50][51] and further emerges in the two-pion contributions via the left-hand cut in γ * γ * → ππ [52][53][54][55][56]. For an isoscalar-photon virtuality q 2 = M 2 ω , M 2 φ , this amplitude is directly related to the three-particle decays of ω and φ, and indeed the KT approach can be generalized to obtain a dispersive representation of the e + e − → 3π cross section, see Sect. 2 for a short review. The ω and φ resonance peaks thus constitute the most conspicuous features of the cross section, but for the HVP integral also the off-peak regions need to be controlled, with QCD determining, via the Wess-Zumino-Witten anomaly [57,58], the normalization in terms of the pion decay constant F π [59][60][61], and, in terms of the KT equations, the ππ rescattering among the final-state pions.
In the 2π channel the average is dominated by experiments using the initial-stateradiation (ISR) technique, while data from energy-scan experiments [62][63][64][65][66][67] are consistent but currently less precise. For the 3π channel this situation is reversed, with only a single ISR data set, which, in addition, only covers the energy region above the φ [68]. Instead, the low-energy region including the ω and φ resonances has been most precisely measured by the Novosibirsk experiments SND [69][70][71][72] and CMD-2 [63,[73][74][75]. For completeness, we will also consider earlier data from DM1 [76], DM2 [77], and ND [78]. The main part of the paper is then devoted to the fit systematics to the various data sets, as detailed in Sect. 3, before working out the consequences for HVP in Sect. 4 and summarizing our findings in Sect. 5.
2 Dispersive representation of the γ * → 3π amplitude Neglecting the mass of the electron, the HVP contribution to (g − 2) µ can be expressed as [22,23] with α = e 2 /(4π), the kernel function 2) and the hadronic cross section 3 R had (s) = 3s 4πα 2 σ(e + e − → hadrons). (2.3) Since higher-order iterations of HVP do become relevant [79,80] (at next-to-leading order, this issue arises, at least in principle, even for HLbL [81]), conventions for the radiative corrections need to be specified. The hadronic cross section is to be understood including final-state radiation (FSR), but with ISR and vacuum polarization (VP) removed ("bare" cross section). This issue of radiative corrections is most severe for the 2π channel, and therein for the ISR data sets, but as demonstrated in [82] the corrections are now known sufficiently accurately that they cannot account for the (g − 2) µ anomaly. For the 3π channel, we have s thr = 9M 2 π in (2.1) and the radiative corrections to the cross section are mainly of conceptual nature. Strictly speaking, a dispersive representation of the γ * → 3π amplitude is only valid in pure QCD, so that, in principle, all photon contributions including FSR should be removed before the fit and only afterwards added again in a perturbative way. In the case of the 2π channel [36], this strategy was indeed carried through in the context of a scalar-QED approximation. For the 3π channel, the full HVP contribution is more than an order of magnitude smaller, so that the total size of the 3πγ final state would be naively estimated at the level 0.3 × 10 −10 , which by itself is borderline relevant at the current level of accuracy. However, since FSR is automatically included in the cross sections provided by experiment, the actual effect only concerns a possible distortion of the fit due to subtracting and adding the FSR contribution, which will be even smaller and therefore neglected here. In contrast, the VP removal does become relevant at the current level of accuracy, mainly because of the resonance enhancement in the vicinity of ρ, ω, and φ, which shifts the pole position, see Sect. 3.5, and modifies the spectral function. When provided by experiment, we use the bare cross section directly, otherwise we apply the VP routine from [25]. To check the sensitivity to this correction, we also constructed an independent VP function based on the 2π fit from [36] as well as the 3π cross section from the present paper, so that major deviations to the full VP only start in the vicinity of the φ, where the KK channels become relevant. We can therefore check the 3π contribution self-consistently up-to-and-including the ω peak, producing a difference of less than 0.1 × 10 −10 in the HVP integral. Accordingly, we conclude that the details of the VP routine lead to a negligible effect as well.
The dispersive representation that we fit to the bare cross section is constructed along the following lines, see [17,18,50] for more details. First, the cross section is given in terms of the γ * → 3π amplitude F(s, t, u; q 2 ) according to with integration boundaries and The amplitude itself is defined by the matrix element of the electromagnetic current j µ with q = p + + p − + p 0 and kinematics The constraints from analyticity and unitarity are most conveniently formulated in terms of the partial-wave amplitudes [83] F(s, t, u; with derivatives of the Legendre polynomials P ′ l (z s ). Since higher partial waves are completely irrelevant below the ρ 3 (1690) resonance [3,51] (see App. B for an estimate of the F -wave contribution), the discontinuity equation reduces to where δ(s) refers to the ππ P -wave phase shift. This is where, in a model-independent way, the information about the ρ(770) enters. The KT equations define, iteratively, the solution of (2.10) in terms of dispersion integrals involving the Omnès function [84] together with deformations of the integration contour necessitated by the decay kinematics. We solve the KT equations with the ππ phase shift recently extracted from the e + e − → π + π − channel [36] and a cutoff parameter Λ 3π = 2.5 GeV. Variations of these input quantities prove irrelevant compared to other sources of systematic uncertainties. For a given q 2 , the KT equations determine the s-dependence of the partial-wave amplitude f 1 (s, q 2 ), but the overall normalization a(q 2 ) is not predicted. At q 2 = 0 it is determined by the low-energy theorem, at q 2 = M 2 ω , M 2 φ it is related to the ω, φ → 3π decay widths, and in general it can be extracted from a fit to the e + e − → 3π cross section. We take essentially the same parameterization as in [17,18] constructed in such a way as to fulfill the low-energy constraint from the chiral anomaly, preserve analyticity of F(s, t, u; q 2 ), and be flexible enough to describe the data up to 1.8 GeV. For the exact relation between the partial wave f 1 (s, q 2 ) and its q 2 -dependent normalization a(q 2 ), see [18,50]. The significance of the individual terms is as follows: the subtraction constant α A is determined by the chiral anomaly (corrected by quark-mass renormalization) [49,85], (2.13) The function A includes resonant contributions via , (2.14) where V = ω, φ, ω ′ (1420), ω ′′ (1650). The energy-dependent widths Γ ω/φ (q 2 ) of the ω/φ mesons include all the main decay channels, in particular, the phase space for the 3π decay channels is calculated including 3π rescattering as well [3] and due to the ω → π 0 γ channel the integration threshold is s thr = M 2 π 0 . For ω and φ, the missing channels account for about 2% of the width, which is remedied by a simple rescaling of the partial widths (in [17,18] the missing ω → π + π − and φ → ηγ were also considered explicitly, leading to virtually identical results). As before, the parameters for ω ′ and ω ′′ are taken from [86], assuming a 100% branching ratio to 3π, but for ω and φ we now allow mass and width to vary: with VP removed, noticeable differences to the PDG emerge, see Sect. 3.5, which is expected since the PDG parameters subsume radiative effects.
Finally, the conformal polynomial in (2.12)  that we interpret as a normalization-type uncertainty and therefore assume to be 100% correlated.
accounts for non-resonant effects. The inelastic threshold s inel is set to 1 GeV 2 motivated by the nearby KK threshold, the second parameter to s 1 = −1 GeV 2 . Further constraints are implemented to remove the S-wave cusp in the polynomial and to ensure that the sum rule is fulfilled exactly. In [17,18] we also introduced further parameters to be able to impose a faster asymptotic behavior of the imaginary part as required for the dispersive description of the pion transition form factor, but since this impaired to some extent the description of the cross section, here, we only consider these additional constraints to estimate systematic uncertainties.
3 Fits to e + e − data

Data sets and unbiased fitting
We start with a brief summary of the data sets that we will include in our analysis, see Table 1. For all data sets the statistical errors are given in diagonal form, with the implica-tion that correlations are negligible at least at the quoted level of uncertainty. In contrast, the treatment of the systematic uncertainties is more ambiguous, since assumptions need to be made on the correlations between data points. Some sources of systematic uncertainty are, by definition, 100% correlated, these are normalization uncertainties for instance due to the luminosity measurement and the detection efficiency, but other systematic effects may well be localized in certain energy regions and therefore should not be considered fully correlated. To follow the experimental documentation as closely as possible, we consider a systematic error of normalization-type origin whenever given as a percentage, otherwise, we treat that uncertainty as a diagonal error. Note that this distinction mainly affects the SND data sets, while for the other energy-scan experiments all systematic errors are given as a percentage. The exception is the ISR data set from BaBar, but [68] states explicitly that the systematic errors for different mass bins are fully correlated. These details are important to monitor a potential bias in the fit. Most importantly, a χ 2 -minimization with an empirical full covariance matrix V(i, j) including a normalization uncertainty, will converge to a solution that is biased towards a lower value than expected due to the fact that smaller data values are assigned smaller normalization uncertainties. This D'Agostini bias was first observed in [87]. It becomes increasingly severe for large normalization uncertainties and/or a large number of data points, so precisely when there is a normalization uncertainty in an experiment that is 100% correlated among all data points. In addition, in a global fit of several experiments a bias that may occur in the combination needs to be avoided. We follow the iterative fit strategy proposed by the NNPDF collaboration [88] to eliminate the bias, which is based on the observation that the normalization uncertainties should be proportional to the true value rather than the measurement. In this manner, the modified iterative covariance matrix is given as where V stat (i, j) is the statistical covariance matrix and the systematic covariance matrix V syst (i, j) is determined by multiplying the normalization factors with the fit function f n (x i ) in each iteration step rather than the data. The empirical covariance matrix can be chosen as the initial guess, with expected rapid convergence to the final solution.
In the fit to the data sets in Table 1 we only encounter either fully correlated or diagonal errors. We follow [88] and treat the uncorrelated systematic errors on the same footing as the statistical ones. For a single experiment one would therefore expect that the central values obtained in a fit with diagonal errors only should be close to the central values of the full fit, otherwise, one would need to understand better the role of the correlations. In the following, we will thus consider both diagonal and full fits to monitor whether significant differences arise.   Table 2: Fits to the combination of SND data sets [69][70][71][72], for diagonal errors and full covariance matrices. p conf denotes the number of free parameters in the conformal polynomial. All errors refer to fit uncertainties only.

Fits to SND
As the first set of fits we consider the SND data sets [69][70][71][72]. The results are summarized in Table 2, both for diagonal errors only and including correlations as described in the previous section. In each case we consider variants of the fits with p conf = 2 . . . 4 free parameters in the conformal polynomial and at this stage display only the fit uncertainties, with systematic uncertainties of the dispersive representation to be added later.
The results in Table 2 show that the main effect of the correlations is an increase in the uncertainty, within the fit statistics the central values agree with the diagonal fit. However, we also note that the description of the data becomes worse, which can be remedied to some extent by increasing p conf . While the diagonal fit proves very stable to variations of p conf , we observe that when including the correlations the central value increases with p conf , balancing the reduction in the central value compared to the diagonal fit in the variant with p conf = 2, the smallest for which a reasonable fit can be obtained.
We note that the treatment of the systematic uncertainties, closely following experiment as specified in Sect. 3.1, is critical to obtain consistent fits. If all systematic uncertainties were assumed to be fully correlated, the fit iteration would not even converge or, when restricted to a subset of the data, lead to a significant downward bias.

Fits to CMD-2 and BaBar
The CMD-2 data sets mainly cover the resonance regions, with [63] scattered around the ω peak and [73][74][75] around the φ. To be able to perform fits to the whole energy region up to 1.8 GeV and thus facilitate the comparison to the SND fits we combine the CMD-2 data with the BaBar data set [68], which starts directly above the φ and covers the remainder. We do not find acceptable fits for the naive combination of all these data sets. To isolate the reason we perform two separate fits, first, to [63,73,74] and BaBar [68], as given in Table 3, as well as [63,75] and BaBar [68], see Table 4. This strategy is motivated by the suspicion that inconsistencies among the CMD-2 data sets arise in the vicinity of the φ, which the separate consideration of the data sets covering this region should be able to corroborate.
In all cases we see that the diagonal and full fits are well compatible, so that the treatment of correlations becomes less of a concern than for the SND fits. However, we find that the fit quality is quite poor: while for the [73,74] φ data sets the fits with p conf = 4 might still be considered acceptable, this is certainly not the case for [75], even though also in this case higher orders in the conformal expansion do yield some improvement of the χ 2 . The most relevant discrepancies in the fit results concern the ω coupling c ω , which is significantly smaller in Table 4 despite being based on the same data set in the ω region, leading to the overall much lower HVP integral, as well as the φ mass.  in Table 3 prefer a value around M φ = 1019.25(4) MeV, while the fits in Table 4 point to M φ = 1019.09(3) MeV, suggesting that inconsistencies in the φ region are compensated elsewhere in the fit, thus the change in c ω . From the mass shifts discussed in App. A, together with the PDG φ mass, we would expect a fit value M φ = 1019.20 MeV, in perfect agreement with Table 2, largely consistent  with Table 3, but clearly at odds with Table 4. Since within uncertainties the ω masses are consistent among the three fits, this suggests as a remedy to include energy-calibration uncertainties in the context of [75], in analogy to the energy rescalings found necessary in the case of the 2π channel [36]. In fact, [75] includes three different scans, and separate fits to each of them reveal that the first two yield φ masses in the expected range, while the third one differs, leading to the lower mass in Table 4. Accordingly, we apply a rescaling to the data of the third scan only. The fit prefers a rescaling around ξ ∼ 10 −4 , well in line with potential uncertainties of the energy calibration. Including ξ as an additional parameter in the fit indeed leads to a mild improvement in the χ 2 . However, removing this tension in M φ by no means renders the resulting fits statistically acceptable. Inspection of the contribution to the χ 2 from each data point shows that a by far disproportionate amount originates from the last few points of each scan of [75] [63,75] and BaBar [68], with modifications to [75] as described in the main text.
which the cross section drops below 5 nb. In the end, we have to conclude that these points cannot be described in a statistically acceptable way with our dispersive representation.
To demonstrate the huge impact on the fit, Table 5 gives the results when these critical points are removed. The fit is clearly still not perfect, but at least comparable in quality to Table 3. Accordingly, we believe that there is reason to suspect some additional systematic uncertainty in the off-peak cross sections from [75] and therefore will only consider the reduced data set as in Table 5 in the following (denoted by CMD-2 ′ ).

Combined fits
Our final preferred fit is shown in Table 6, including all data sets listed in Table 1 except for the DM2 data [77], which disagree with both the BaBar [68] and the SND [72] data especially in the vicinity of the ω ′′ (1650). We also considered fits dropping the CMD-2 data set [75] altogether, see Table 7, but the overall effect is relatively minor, depending on the fit variant at most 0.2 × 10 −10 in the final (g − 2) µ integral. In all cases the χ 2 is significantly worse than in the separate fits discussed in the previous sections. Since there are now several experiments covering the same energy region, this is an indication of the degree of consistency among the various data sets. To account for these   (9) 10 4 × ξ 1.9(7) 1.8 (7) 1.8 (7 [76], and ND [78]. inconsistencies we follow the PDG prescription [86] and inflate the fit errors by the scale factor S = χ 2 /dof, (3.4) which increases uncertainties by about 20% compared to the fit errors given in Table 6. The systematic uncertainties are dominated by the degree of the conformal polynomial, while the uncertainties from ππ phase shifts and cutoff parameters are negligible in comparison. We adopt the full results for p conf = 3 as our central value (as this gives the best fit), but keep the maximum differences to p conf = 2, 4 as a source of systematic uncertainty. In addition, we perform fits in which the imaginary part of the conformal polynomial in (2.15) is constrained to behave as q −3 asymptotically, and include the observed variation as another source of systematics, see Table 8 for this last set of fits. As alluded to earlier, the fit quality deteriorates when imposing this additional constraint on the conformal polynomial, and therefore the full variation over all fit variants would likely be an overestimate of the systematic uncertainty. To gauge the impact, we take the average change for p conf = 2, 3, 4 separately, and add the result in quadrature to the systematic error from the variation in p conf . 4 The final fit is illustrated in Fig. 1

Extracting ω and φ masses
Our final result for the ω and φ parameters is one needs to keep in mind that these parameters subsume radiative effects, where the expected corrections are worked out in App. A. For the φ, the expectation is that the fit of large shifts in the fit parameters compared to p conf = 3 and sizable cancellations among the terms in the conformal polynomial. We still include this fit in the estimate of the systematic uncertainties, otherwise, the systematic errors of the final results given in Sects. 3 Table 6, but with Im C p (q 2 ) ∼ q −3 asymptotically. In view of the expected downward shift of 0.13 MeV, our analysis thus supports the 3π number from [71], while the agreement with the PDG average is entirely coincidental. As argued in [36], the π 0 γ value is likely affected by an unphysical phase in the extraction, but our analysis shows that a similar effect does not occur in 3π. Therefore, our analysis compounds the tension with the VP-subtracted ω mass as extracted from the 2π channel, M ω = 781.68 (9) Figure 1: Fit to the e + e − → 3π data sets as listed in Table 1 (with VP removed everywhere). The black band includes the fit uncertainties only, while the gray band represents the total uncertainty, including the systematics of the dispersive representation. For most energies the two uncertainties are of similar size, so that the difference is hardly visible on the logarithmic scale.
where the systematic errors are estimated as in Sect. 3.4. As a cross check we have also performed a fit to the data combination of [25] instead of the data directly, leading to almost the same central value with slightly larger uncertainties. The latter is likely related to the fact that although, as expected, we had to remove two bins (# 49 and # 52) corresponding to (nearly) vanishing cross sections to have the fit iteration converge, no further changes were applied to the combination, so that some of the potentially problematic points we identified in [75] could still impact the fit. The final result (4.1) agrees well with a 3π µ | ≤1.8 GeV = 46.2(1.5) × 10 −10 [24], besides a corroboration of the central value the QCD constraints also allow for a reduction of the uncertainty. The difference to a 3π µ | ≤1.8 GeV = 47.7(9) × 10 −10 [25] is mainly due to the interpolation applied to the data. We reproduce the central value with a linear interpolation of the bins of [25], while higher-order interpolations, as well as the dispersive fit (4.2), move the central value towards (4.1). 5 Our analysis does not support values as low as a 3π µ | ≤2.0 GeV = 44.3(1.5) × 10 −10 [91], which is based on a Breit-Wigner description of ω and φ. Finally, we remark that for the threshold region we find a value a 3π µ | ≤0.66 GeV = 0.019 × 10 −10 nearly twice as large as the estimate from [92] based upon a combination of the Wess-Zumino-Witten action and vector meson dominance [93,94]. Indeed, it was observed in [92] that this model underestimates the lowest-energy data points.
In combination with the 2π channel from [36] we obtain for the HVP contribution that has been evaluated imposing analyticity and unitarity constraints Even assuming for the HLbL contribution a value as conservative as a HLbL µ = 10(4) × 10 −10 , this result thus reaffirms the (g − 2) µ anomaly at the level of 3.4σ. 6

Summary
We have presented a detailed analysis of the 3π contribution to HVP, including constraints from analyticity and unitarity as well as the low-energy theorem for the γ * → 3π amplitude. Similarly to the 2π analysis of the pion vector form factor [36], the main motivations are, first, to see if a global fit subject to these constraints reveals inconsistencies in the data, and, second, derive the corresponding error estimate for the contribution to (g − 2) µ . Given that this method is complementary to a direct integration of the data, where potential inconsistencies are addressed by a local error inflation, such global fits that incorporate general QCD constraints should increase the robustness of the SM prediction.
We find that most data sets can be fit satisfactorily with our dispersive representation, the exception being several points above the φ resonance from [75]. Fortunately, the impact on the final HVP integral is minimal, due to the suppression of the cross section in this region, which could also enhance the relative importance of systematic effects in the data. Otherwise, in the 3π channel there is no tension between two high-statistics data sets, such as BaBar and KLOE in the 2π case, but the scale factor of the global fit, indicating overall consistency of the data base, is actually larger than in 2π. In addition, the main contribution, from the 3π cross section in the vicinity of the ω, is dominated by a single experiment [71]. For these reasons, a new high-statistics low-energy measurement in the 3π channel would be a highly welcome addition to the data base.
The central outcome of our study is (4.1) Together with the 2π channel from [36], the two most important low-energy channels have now been scrutinized including analyticity and unitarity constraints, covering nearly 80% of the HVP integral. Depending on the assumptions for HLbL scattering, the current tension is thus confirmed at the known level around 3.5σ. To make further progress, especially in view of the ultimate precision expected at the Fermilab (g − 2) µ experiment, new data input in particular in the 2π channel is critical. Finally, our analysis exacerbates a tension emerging between the 2π and 3π channels, that is, the extraction of the ω mass. In the 2π channel the ω only contributes via an isospin-violating effect, ρ-ω mixing, but due to the increased statistics the sensitivity is not much below that of the 3π channel. Yet, the ω mass extracted from the 2π channel is substantially lower than the one extracted from 3π. Currently, we are aware of neither a systematic effect in experiment nor an issue with the theoretical extraction that could resolve the tension. Besides improving the HVP contribution to (g − 2) µ , new data could shed light on this puzzle as well.

A Electromagnetic mass shifts
The separation of VP from the full cross section affects the ω and φ pole parameters because the VP function itself involves the corresponding poles, only suppressed by e 2 . The size of the expected shifts can be analyzed analytically in a Bethe-Salpeter multi-channel approach [101,102]. For instance, the ω contribution to the VP function Π(s) becomes Π ω (s) = e 2 s g 2 where the coupling is related to the two-electron width Γ ω→e + e − = e 4 M ω /(12πg 2 ωγ ), i.e. g ωγ = 16.7(2) [86]. Expanding around the shifted pole parameters, one finds the relation M ω is thus expected to be about 0.13 MeV lower than in PDG conventions, while the effect on the width due to Π ω is negligible. The same argument for the φ produces a mass shift Finally, for the ω width there is an additional effect due to ρ-ω mixing, i.e., a higherorder effect enhanced by the small mass difference between ω and ρ. In a vector-mesondominance approximation for the ρ we find the relation with mixing parameter ǫ ω ∼ 2 × 10 −3 [36], and by comparing fits with and without VP we verified that this indeed describes well the observed shift in the ω width.

B Estimate of the F -wave contribution
For q 2 = 0 [51] and q 2 = M 2 ω [3] the impact of F -waves on the γ * → 3π amplitude was shown to be completely negligible below the ρ 3 (1690) resonance, but since we consider virtualities up to q 2 = 1.8 GeV one may ask the question whether the impact of these resonant F -waves can still be ignored. There is little phenomenological information on the ρ 3 πγ * coupling besides the ρ 3 → πω branching ratio. However, the fact that the corresponding ω-dominance estimate from [51] is in line with preliminary results from COMPASS [103] suggests that at least within [0, M 2 ω ] the q 2 -dependence should be approximately described by a(q 2 ). Here, we estimate a potential F -wave contribution by assuming that this approximation remains meaningful up to q 2 = 1.8 GeV.
The decomposition of the amplitude including F -waves becomes [3] F(s, t, u; q 2 ) = F(s; q 2 ) + F(t; q 2 ) + F(u; q 2 ) + P ′ 3 (z s )G(s; q 2 ) + P ′ 3 (z t )G(t; q 2 ) + P ′ 3 (z u )G(u; q 2 ), (B.1) where the scattering angles follow by permuting the Mandelstam variables in (2.8) accordingly. To estimate the ρ 3 contribution, we first establish the connection to a narrowresonance approximation of the P -wave with M 2 ρ → M 2 ρ − iM ρ Γ ρ in the decay region. We can then estimate the ρ 3 contribution as where the coupling constants are set to the values from [51]. Numerically, we find that the interference between P -and F -waves gives a correction around 1% at q 2 = 1.8 GeV, while the pure F -wave contribution is suppressed by another two orders of magnitude. These results confirm the expectation that the ρ 3 (1690) should not become relevant until well above the threshold M ρ 3 + M π ∼ 1.83 GeV where the decay becomes possible.