Two-pion contribution to hadronic vacuum polarization

We present a detailed analysis of e+e− → π+π− data up to s=1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sqrt{s}=1 $$\end{document} GeV in the framework of dispersion relations. Starting from a family of ππ P-wave phase shifts, as derived from a previous Roy-equation analysis of ππ scattering, we write down an extended Omnès representation of the pion vector form factor in terms of a few free parameters and study to which extent the modern high-statistics data sets can be described by the resulting fit function that follows from general principles of QCD. We find that statistically acceptable fits do become possible as soon as potential uncertainties in the energy calibration are taken into account, providing a strong cross check on the internal consistency of the data sets, but preferring a mass of the ω meson significantly lower than the current PDG average. In addition to a complete treatment of statistical and systematic errors propagated from the data, we perform a comprehensive analysis of the systematic errors in the dispersive representation and derive the consequences for the two-pion contribution to hadronic vacuum polarization. In a global fit to both time- and space-like data sets we find aμππ|≤ 1 GeV = 495.0(1.5)(2.1) × 10− 10 and aμππ|≤ 0.63 GeV = 132.8(0.4)(1.0) × 10− 10. While the constraints are thus most stringent for low energies, we obtain uncertainty estimates throughout the whole energy range that should prove valuable in corroborating the corresponding contribution to the anomalous magnetic moment of the muon. As side products, we obtain improved constraints on the ππ P-wave, valuable input for future global analyses of low-energy ππ scattering, as well as a determination of the pion charge radius, 〈rπ2〉 = 0.429(1)(4) fm2.


Introduction
ππ scattering is one of the simplest hadronic reactions that displays many key features of low-energy QCD [1], most prominently approximate chiral symmetry, its spontaneous breaking, and the explicit breaking due to finite up-and down-quark masses. Accordingly, chiral symmetry severely constrains the low-energy scattering amplitude, which can be systematically analyzed in Chiral Perturbation Theory (ChPT) [2][3][4][5] and has been worked out up to two-loop order [6]. However, ππ scattering is not only unique because of its strong relation to chiral symmetry, but in addition exhibits further remarkable properties JHEP02(2019)006 that extend beyond the low-energy region where the chiral expansion applies. Most notably, this includes the fact that the process is fully crossing symmetric and that the unitarity relation, up to center-of-mass energies of nearly √ s = 1 GeV, is totally dominated again by ππ scattering. The resulting constraints from analyticity, unitarity, and crossing symmetry were first formulated systematically in the framework of Roy equations [7] and subsequently analyzed in great detail, ultimately leading to a very precise representation of the ππ phase shifts up to roughly 1 GeV [8][9][10]. Both the methods used in determining these phase shifts and the actual results have had a profound impact on countless more complicated hadronic reactions and decays, such as πK scattering [11,12], πN scattering [13,14], η → 3π [15][16][17], η → ηππ [18], or K 4 decays [19]. However, arguably the most immediate application concerns pion form factors and here especially the vector form factor F V π , given that, in marked contrast to the scalar form factor [20][21][22][23], the onset of inelastic corrections is relatively smooth.
Since the early determinations [32][33][34][35] the experimental situation in e + e − → π + π − has improved considerably [63][64][65][66][67][68][69][70][71][72][73][74][75], but at the same time the required precision of the HVP contribution to a µ has increased further, in particular in view of the anticipated improve-JHEP02(2019)006 ment of the experimental measurement of a µ by a factor 4 at the Fermilab experiment. In this way, a proper treatment of experimental errors and correlations is becoming absolutely critical. This includes radiative corrections, which need to be taken into account properly in order to ensure a consistent counting of higher-order HVP iterations [76,77] (in principle, the same issue arises in HLbL scattering as well [78]). Most current HVP compilations are based on a direct integration of the experimental data [79][80][81] (see [82] for an alternative approach using the hidden-local-symmetry model), wherein conflicting data sets are treated by a local χ 2 inflation. 1 The most consequential such tensions currently affect the BaBar [70,72] and KLOE [69,71,73,75] data sets for the ππ channel, and different methods for their combination then give rise to the single largest difference between the HVP compilations of [80] and [81].
In this paper, we return to the description of the ππ contribution to HVP based on a dispersive representation of the VFF. We first clarify the role of radiative corrections, in particular vacuum polarization (VP), and then derive a global fit function that the form factor needs to follow to avoid conflicts with unitarity and analyticity. In addition to two free parameters in the ππ P -wave, this Omnès-type dispersion relation involves one parameter to account for ρ-ω mixing (plus the ω mass) and at least one additional parameter to describe inelastic corrections in a conformal expansion. First, we study to which extent the resulting representation can be fit to the modern high-statistics data sets, by using an unbiased fit strategy and including the full experimental covariance matrices where available, to provide a strong check of the internal consistency of each data set. As a second step, we address the systematic uncertainties in the dispersive representation and derive the HVP results for various energy intervals. Finally, we provide the resulting ππ P -wave phase shift and the pion charge radius that follow after determining the free parameters from the fit to the e + e − → π + π − cross section data.

Dispersive representation of the pion vector form factor
In this section we review the formalism for a dispersive representation of the pion VFF at the level required for the interpretation of the modern high-statistics e + e − → π + π − data sets. This includes the definition of the pion VFF in QCD, the relation to HVP, and conventions regarding radiative corrections, see section 2.1; the actual dispersive representation including the description of the most important isospin-breaking effect from ρ-ω mixing as well as a term accounting for inelastic contributions, see section 2.2; and a constraint on the size of these inelastic contributions, the Eidelman-Lukaszuk bound, see section 2.3.

Hadronic vacuum polarization and radiative corrections
Hadronic contributions to the anomalous magnetic moment of the muon first arise at O(α 2 ) in the expansion in the electromagnetic coupling α = e 2 /(4π). The leading topology is HVP, shown in figure 1, with hadronic input encoded in the QCD two-point function of JHEP02(2019)006 electromagnetic currents Π µν = ie 2 d 4 xe iq·x 0|T {j µ em (x)j ν em (0)}|0 = (q 2 g µν − q µ q ν )Π(q 2 ), (2.1) where the Lorentz decomposition follows from gauge invariance, the current is defined by 2 j µ em :=qQγ µ q, q = (u, d, s) T , Q = diag , (2.4) where in pure QCD the integral starts at the two-pion threshold, s thr = 4M 2 π . Unitarity relates the imaginary part of Π(s) to the total hadronic e + e − cross section σ(e + e − → hadrons) = α s 4π σ e (s) 1 + 2m 2 e s Im Π(s), (2.5) where σ (s) = 1 − 4m 2 /s. The HVP contribution to the anomalous magnetic moment of the muon can then be written as [89,90] a µ = αm µ 3π 2 ∞ s thr dsK (s) s 2 R had (s), (2.6) where the kernel function iŝ and R had is related to the hadronic cross section by R had (s) = 3s 4πα 2 sσ e (s) s + 2m 2 e σ(e + e − → hadrons). (2.8) We stress that the usual R ratio, defined as the ratio of hadronic to muonic e + e − cross sections, is not what enters the dispersive representation of the HVP contribution: our representation (2.6), with (2.7) and (2.8) as input, is exact. R had (s) and R(s) coincide for 2 As usual in the context of g − 2, we do not include e in the definition of the current. However, we keep Π µν and Π as quantities of O(e 2 ) by including the explicit factor of e 2 in (2.1). 3 For simplicity, in the following we drop the superscript "ren".

JHEP02(2019)006
a tree-level muonic cross section and in the limit s m 2 µ , where of course also the electron mass does not play any role, but for clarity, we have provided above the expression of the HVP contribution without any approximations.
The contribution of the two-pion intermediate state can be expressed in terms of the pion VFF π ± (p )|j µ (2.10) As becomes apparent already from (2.3), for a consistent counting of higher orders in α it is critical that radiative corrections be properly taken into account, otherwise this would induce corrections at the same order as HVP iterations or HLbL scattering. The prevalent convention is that the leading-order HVP include, in the sum over intermediate states in the unitarity relation, not only the hadronic channels but the (one-)photon-inclusive ones. In particular, the lowest-lying intermediate state is no longer the two-pion state, but the π 0 γ state, i.e. s thr = M 2 π 0 . In this way, the HVP input corresponds to infrared-finite photoninclusive cross sections including both real and virtual corrections, but to avoid double counting in next-to-leading-order iterations and beyond each contribution is required to be one-particle-irreducible. This convention has important consequences for the definition of the pion VFF and the corresponding e + e − → π + π − cross section. That is, the cross section to be inserted in (2.8) has to be inclusive of final-state radiation (FSR), but both VP and initial-state-radiation (ISR) effects need to be subtracted. This defines the bare cross section where the running of α, see (2.3), is determined by the full renormalized VP function in the SM, e.g. including the lepton-loop contribution While by means of the above equations the subtraction of VP effects may be taken into account afterwards, the correction of ISR and ISR/FSR interference effects is performed with Monte Carlo generators in the context of each experiment [91][92][93][94]. Accordingly, the two-pion contribution should be understood as the photon-inclusive two-pion channel. This, however, is not directly compatible with our goal to treat the pion VFF dispersively, because this is usually done in the isospin limit, i.e. with photon emission switched off. In order to be able to apply our dispersive treatment of the VFF, we therefore need to extract from the data on the photon-inclusive process the cross section σ(e + e − → π + π − ) in the isospin limit, i.e. with m u = m d and α = 0. While taking this JHEP02(2019)006 limit is unproblematic at first sight, subtleties arise once one realizes that experiments exist only in our isospin-broken world and that any input quantity has to be taken and defined away from the isospin limit. Phrased differently, the actual question one faces is whether it is possible to establish a procedure to extract from the measured photoninclusive cross section σ(e + e − → π + π − (γ)) the one in the isospin limit, where the VFF in pure QCD appears.
A similar question shows up also in other contexts. One case that has been discussed in detail in the literature is the problem of the extraction of the purely strong pion decay constant from the measurement of the decay rate Γ(π → µν µ (γ)). Strictly speaking, an unambiguous and uniquely defined extraction is not possible: any practical and operative definition of a purely strong decay constant extracted from experiment is necessarily convention dependent. More precisely, it depends on how one defines the strong isospin limit and on a matching scale, as explained in detail in [95,96] from the perspective of an effective-field-theory, perturbative approach, and in [97] from a non-perturbative point of view. For the pion decay constant it has been shown [96] that the scale dependence is very weak, mainly thanks to the smallness of α and the logarithmic dependence on the scale, and that the extraction of the pion decay constant from experiment is indeed accurate at the claimed accuracy, of course barring the choice of absurd values of the matching scale.
An example where the scale dependence in defining a purely strong quantity cannot be neglected concerns proton-proton scattering. Here, the scale-independent photoninclusive scattering length a C pp = −7.8063(26) fm [98] differs significantly from the scaledependent photon-subtracted one [99,100], depending on the choice of scale e.g. a pp = −17.3(4) fm [101]. In this case, the size of the effect is enhanced by the interference of the Coulomb interaction with the short-distance part of the nuclear force, and virtual photons could only be subtracted consistently everywhere, including the running of operators, if the underlying theory were known [96,102]. This situation should be contrasted with perturbative systems, e.g. the extraction of the pion-nucleon scattering lengths from pionic atoms [103][104][105], where in principle the same ambiguities related to the removal of QED effects appear, but, without such an enhancement mechanism, the resulting scale dependence can be neglected at the level of the experimental accuracy.
For the case of the e + e − → π + π − (γ) cross section the situation is completely analogous to that of F π : in principle, the purely strong VFF cannot be extracted in an unambiguous way from data, but one may hope that a convention-dependent extraction (and corresponding definition) of such a strong VFF only shows a very weak dependence on this arbitrariness and can be taken as a good approximation to a purely strong VFF. The problem has been analyzed in the literature mainly with the help of scalar QED and extensions thereof that include resonance exchanges [93,94,[106][107][108][109][110][111]. In these models there is no ambiguity in the extraction of the cross section σ(e + e − → π + π − ), but this happens at the price of losing model independence. In either case, these studies indicate that at the present level of accuracy scalar QED describes reasonably well the behavior of the observed FSR: the relation established within this model between the cross section without photon emission and the fully inclusive one is likely to be sufficiently accurate for our purposes,
is absolutely essential for our purposes: otherwise the dispersive constraints to be discussed in the next section would not apply.
There are other issues related to radiative corrections which have also been discussed in the literature, in particular ρ-ω and ρ-γ mixing [112], in the context of a Bethe-Salpeter approach for the coupled-channel system of e + e − , π + π − , and 3π, see [113,114]. The main result is that additional corrections from ρ-γ mixing [112] only become relevant if an attempt is made to define external ρ states, as required for estimates of isospin-breaking corrections in the interpretation of τ → ππν τ data to be used as input for HVP [115,116], but in the e + e − → π + π − channel full consistency is ensured as long as the same pure-QCD form factor F V π that determines the bare cross section (2.14) defines, self-consistently, the π + π − contribution to the VP function Π SM in its extraction from experiment (2.11). In practice, we find that the VP routines applied in the modern experiments are sufficiently close to such a fully self-consistent solution that we can use the bare cross sections as provided by experiment. 5 Accordingly, the physical FSR-inclusive cross section takes the form where the VP function has been expressed in terms of the running coupling α(s). Unfortunately, the common procedure in the literature amounts to absorbing a factor α(s)/α into the definition of the form factor, see [63][64][65][66][67][68][69][70][71][72][73][74][75], but in these conventions we could not 4 A similar factorization assumption has recently been used in order to extract the η → 3π differential decay rate in the isospin limit from the measured one [17]. 5 In fact, if the normalization is determined from e + e − → µ + µ − , the resulting cross section is automatically bare because VP drops out in the ratio. This applies to the BaBar [70,72] and KLOE12 [73] data sets.

JHEP02(2019)006
formulate the dispersive constraints. For this reason, we do not use the results for F V π (s) provided by each e + e − experiment, but rather the bare cross section in order to reconstruct the actual QCD form factor.

Omnès representation of the form factor
In the following, we present the dispersive representation of the pion VFF F V π (s) as put forward in [33,34]. In particular, we treat the form factor in pure QCD and include the most important strong isospin-breaking effect from the mixing into the 3π channel. In the isospin limit, F V π (s) is an analytic function of s, apart from a branch cut in the complex s-plane that lies on the real axis, s ∈ [4M 2 π , ∞), and is dictated by unitarity. The form factor is real on the real axis below the branch point 4M 2 π , hence it fulfills the Schwarz reflection principle. We parametrize the pion VFF as a product of three functions, where is the usual Omnès function [117] with δ 1 1 (s) the isospin I = 1 elastic ππ phase shift in the isospin-symmetric limit. The Omnès function alone is the solution for the VFF in the isospin limit and disregarding inelastic contributions to the unitarity relation. Therefore, the quotient function F V π (s)/Ω 1 1 (s) is analytic in the complex s-plane apart from a cut on the real axis starting at s = 9M 2 π . The factor G ω accounts for ρ-ω mixing, the most important isospin-breaking effect, which becomes enhanced by the small mass difference between the ρ and ω resonances. The full parametrization implements the correct threshold behavior of the discontinuity, i.e. the right-hand cut starting at 9M 2 π opens with the fourth power of the center-of-mass momentum [33]. In practice, it would even be possible to replace G ω (s) by g ω (s) with almost no observable effect in the energy range of interest, in particular, due to the strong localization around the ω resonance, the imaginary part of g ω (s) below threshold is tiny. We still use the dispersively-improved variant (2.18) to remove this unphysical imaginary part altogether and have the threshold behavior correct, but stress that if the difference became relevant, this form would not suffice to go beyond the narrow-width approximation. For that also the spectral shape would need to be improved [53,56]. For completeness, we remark that JHEP02(2019)006 while in general P -wave phase space predicts a behavior proportional to (s−n 2 M 2 π ) 3(n−1)/2 , the leading term vanishes for n = 3, giving rise to the extra power in (2.18).
The remaining function G N in (s) is analytic in the complex s-plane with a cut on the real axis starting at s = 16M 2 π . It takes into account all further inelastic contributions to the unitarity relation. We describe it by a conformal polynomial (2.20) where the conformal variable is and we consider inelasticities only above s in = (M π 0 + M ω ) 2 , since 4π inelasticities are extremely weak below, see section 2.3. The conformal polynomial generates a branch-cut singularity at s = s in and in the variant (2.20) does not modify the charge, G N in (0) = 1. We also require the cut to reproduce P -wave behavior at the inelastic threshold, i.e. close to s in the function G N in (s) has to behave like (s in − s) 3/2 Hence, in order to have a vanishing coefficient of the √ s in − s term, we impose In summary, our parametrization of the form factor fulfills all requirements of analyticity and unitarity, including explicitly the 2π and 3π channels and inelastic corrections in a conformal polynomial with threshold dictated by phenomenology. We expect this representation to be accurate as long as the conformal polynomial provides an efficient description of inelastic effects, conservatively estimated below √ s = 1 GeV. As main input, we require the elastic ππ P -wave phase shift δ 1 1 (s), see section 3, while the isospin-breaking and inelastic corrections are parametrized in terms of the ω parameters ( ω , M ω , and Γ ω ) and c k and s c in the conformal polynomial, respectively.

Inelastic contributions and Eidelman-Lukaszuk bound
In [118,119], a generalization of Watson's theorem [120] was derived that amounts to a constraint on the difference between the phase of the VFF and the elastic ππ scattering phase shift, the Eidelman-Lukaszuk (E L) bound: E L bound with ι 1 = 0.05(4) Figure 2. E L bound on the difference between the phases of the pion VFF and the elastic ππ P -wave. The bound uses the data compilation of [118] for the cross section ratio r, and an elasticity parameter calculated with ι 1 = 0.05 (4). The smaller black error bars indicate the uncertainty due to r, the larger gray error bars the uncertainty due to ι 1 .
where ψ denotes the phase of the form factor, F V π (s) = |F V π (s)|e iψ(s) , η 1 is the ππ elasticity parameter, defined by the expression for the ππ P -wave amplitude and r is the ratio of non-2π to 2π hadronic cross sections r = σ I=1 (e + e − → hadrons) σ(e + e − → π + π − ) − 1 (2.26) in the isospin I = 1 channel. For r < 1, the bound (2.24) implies η 1 ≥ (1 − r)/(1 + r), resulting in [119] sin 2 (ψ − δ 1 1 ) ≤ With a given input for the elasticity parameter η 1 , the bound (2.24) usually provides a much stronger constraint than (2.27), but the latter shows that a non-trivial bound arises as soon as r > 0 irrespective of η 1 . We use a representation of the elasticity parameter from the ππ Roy-equation analysis [8,10] with s a = (1 GeV) 2 and ι 1 = 0.05 (4). With the experimental input on r from [118], we obtain the bound on the phase difference shown in figure 2. Using the parametrization (2.28) for small values of ι 1 , the bound can be conveniently written as

JHEP02(2019)006
where the negligible O(ι 2 1 ) term is negative as long as r 2 < 3. The E L bound provides an important constraint on the parameters c k of the conformal polynomial that we use to describe the inelastic contributions.
We note that in contrast to the value of ι 1 = 0.05(5) from [8,10], we vary the parameter in a slightly smaller range in order to exclude a vanishing value of ι 1 , which would correspond to η 1 = 1, while a non-zero value of r always implies η 1 < 1. 6 In principle, very small values of ι 1 could be excluded by considering particular channels, such as the π 0 ω intermediate state [50] that motivates the functional form (2.28). At slightly higher energy, theKK channel opens, which gives a rather small contribution to the inelasticity [11,50], but also shows that at some point η 1 = 1 is excluded by data. Here, we motivate the lower bound on ι 1 directly through the fits to the e + e − data: if ι 1 is chosen too small, the conformal polynomial becomes constrained too much, resulting in an unacceptable fit quality. In this way, the e + e − data themselves imply that the inelasticity cannot be arbitrarily small. Our range covers those values for which the fits are still acceptable.

Input for the phase shift
The central input in our representation of the pion VFF is the elastic P -wave ππ scattering phase shift. We use the solution of the Roy-equation analysis [8,10]. The parametrization of the phase shift of [10] depends on 27 parameters, but most of them concern elasticity parameters or input from Regge theory for the asymptotic region, both in the P -wave and the other amplitudes related by crossing symmetry. The (elastic) P -wave phase shift itself, below 1.15 GeV, only involves two free parameters, which can be identified with its values at s 0 = (0.8 GeV) 2 and at s 1 = (1.15 GeV) 2 , whose current estimates are [10]  15 GeV) 2 as adopted in [10] the mathematical properties of the equations dictate that there be exactly two free parameters [121,122]. In our description of the VFF, the values of the phase at s 0 and s 1 enter as fit parameters, while the values (3.1), derived from previous analyses of the VFF, only serve for comparison and as starting values in the fit. All the remaining 25 parameters of the Roy solution will be varied within the ranges given in [10] and treated as a source of systematic uncertainties in our description. The central solution for the phase shift is shown in figure 3 together with an uncertainty band generated by varying these 25 parameters. At energies above 1.3 GeV, the ππ phase shift is not as well known as in the lowenergy region shown in figure 3. However, this uncertainty will not have a strong impact on the low-energy description of the form factor. We estimate this uncertainty by studying different prescriptions for the high-energy continuation of the phase shift. Asymptotically, Elastic ππ P -wave phase shift δ 1 1 Figure 3. Elastic P -wave ππ scattering phase shift δ 1 1 from the solution of the Roy equations [10]. We show the central phase solution with an uncertainty band generated by varying all parameters apart from the phase values at s 0 and s 1 .
we assume that the phase shift reaches [33,[123][124][125][126][127][128][129] so that the Omnès function behaves asymptotically as Ω 1 1 (s) s −1 . For our central phase solution, we use the simple prescription [130] with s a = (1.3 GeV) 2 , and we compare to the prescription Continuations of the ππ P -wave phase shift δ 1 From the ππ scattering phase shift δ 1 1 we calculate the Omnès function (2.17), with squared absolute value shown in figure 5. The uncertainty band is again generated by varying all the parameters of the Roy solution apart from the phase values at s 0 and s 1 , which for this plot we have fixed at the central values (3.1). Note that although the Omnès function already closely resembles the pion VFF, the uncertainty of |Ω 1 1 | 2 shown in the plot will not translate directly to |F V π | 2 , because for the description of the VFF δ 1 1 (s 0 ) and δ 1 1 (s 1 ) will not be fixed but enter as fit parameters. The fact that only two free parameters are allowed in the description of the phase shift emphasizes the stringent constraints that follow from ππ Roy equations, ensuring in each step that the solution for the phase shift is consistent with analyticity, unitarity, and crossing symmetry. This is the crucial advantage over using a phenomenological parameterization of the phase shift instead, based on which a confrontation of the VFF data with these general QCD properties would not be possible.

Fits to e + e − data
In this section, we first describe in section 4.1 the parameters in our representation of the pion VFF. They are either fit to experimental data or treated as sources of systematic uncertainties in the theoretical description. In section 4.2, we give an overview of the available data sets and describe the procedure that we use to avoid bias in the fit. In Omnès function from elastic ππ P -wave phase shift perform fits to combinations of the data sets. In section 4.5, we compare the fit result for the ω mass with extractions from other channels and discuss the observed tension.

Fit parameters and systematic uncertainties
The representation of the pion VFF (2.16) is given by a product of three functions. Each of them contains parameters that we fit to data from e + e − → π + π − experiments. First, the Omnès function Ω 1 1 contains two free parameters, the values of the elastic ππ scattering P -wave phase shift at two points, δ 1 1 (s 0 ) and δ 1 1 (s 1 ). Second, the function G ω involves the ρ-ω mixing parameter ω as a free parameter, while the ω mass will either be taken as an input or considered a fit parameter as well. Third, the function G N in describing inelastic contributions contains the N − 1 fit parameters c k of the conformal polynomial (c 1 is constrained by (2.23)). Finally, we will also consider fit variants in which we allow for an experimental uncertainty in the energy calibration, which we will implement by rescaling the energies of the data points constrained by the expected calibration uncertainty of each experiment in the vicinity of the ρ peak. For a single experiment, this strategy produces a similar effect as fitting the ω mass, but for the combined fits it allows us to separate a single ω mass as determined from the e + e − → π + π − fits from variations among the different experiments that might be attributed to the energy calibration.
All other parameters in the form factor representation are treated as sources of systematic uncertainties. First, the 25 additional parameters in the solution of the Roy equations for the phase shift δ 1 1 are varied independently within the ranges estimated in [10], with the JHEP02(2019)006 exception of ι 1 , which determines the elasticity (2.28) and does not only appear in the phase shift but also in the E L bound (2.24). This parameter is varied within ι 1 ∈ [0.01, 0.09], as explained in section 2.3. A second source of systematic uncertainty concerns the continuation of the phase shift to energies above the validity of the Roy equations as described in section 3. If not fit to e + e − → π + π − data, the omega mass is taken as an input from the PDG [132] M ω = 782.65 (12) MeV. (4.1) Since we do not observe any improvement of the fits by letting the ω width float, we keep it as an input [132] Γ ω = 8.49 (8) MeV.
Next, in the conformal polynomial, the point s c that is mapped to the origin is a free parameter. It should be taken sufficiently far from any branch cuts. We vary it in the range and treat it as another source of systematic uncertainty. Finally, the order N of the conformal polynomial is varied between N = 2 and N = 6.

Data sets and unbiased fitting
In our fits of the pion VFF, we take into account the high-statistics time-like data sets from the e + e − experiments. On the one hand, there are the results from the energy-scan e + e − → π + π − experiments SND [65,66] and CMD-2 [63,64,67,68] at the VEPP-2M collider in Novosibirsk. On the other hand, the so-called radiative return measurements run at a fixed e + e − energy and vary the π + π − energy by making use of ISR in the process e + e − → π + π − γ. These experiments are BaBar [70,72] at SLAC, KLOE [69,71,73,75] at the Frascati φ-factory DAΦNE, and BESIII [74] at the BEPCII collider in Beijing.
In addition to these time-like data sets, there is also some experimental information on the space-like form factor available from the scattering of pions off an electron target, performed by the F2 experiment [133] at Fermilab and by NA7 [134] at CERN. Although we have checked consistency of the fit with the extraction of the space-like form factor from e − p → e − π + n data by the JLab F π collaboration [135][136][137][138], we do not use these data in our fits because of their remaining model dependence due to the extrapolation to the pion pole.
For all the data sets that we use in the fits, the experimental uncertainties are split into statistical and systematic errors, see table 1 for an overview. In the case of the space-like data sets and the energy-scan e + e − experiments, the statistical uncertainties are assumed to be uncorrelated between the data points. The systematic errors in general are multiplicative uncertainties similar to overall normalization errors. If fits to data with this type of uncertainties are performed by minimizing a χ 2 function that is constructed with the naive covariance matrix one usually introduces a bias, as first observed by D'Agostini [140]. The bias can be severe especially when combining different data sets with normalization uncertainties. We use an iterative method to avoid this bias as proposed by the NNPDF collaboration [141]. To this end, we define a systematic covariance matrix for relative values and use the following covariance matrix in the χ 2 function: i.e. the relative systematic covariance is weighted by the values of the fit model and not the data. We assume some initial value for the model parameters and iterate the fit with a new covariance matrix constructed using the model function of the previous fit iteration. The iterative fit converges rapidly, typically after only a couple of iteration steps.
In the case of the space-like data sets, our fit function f (x i ) is the squared modulus of the form factor |F V π (s i )| 2 at the center-of-mass squared energies s i of the data points. For all time-like data sets, we use the bare cross sections, which are already undressed of VP effects, and we correct for FSR effects as explained in section 2.1 to relate the bare cross section to the form factor in pure QCD. In contrast, the VFF data directly provided by experiment still contain VP effects and therefore cannot be consistently fit with our QCD-only parametrization. Hence, in the case of the energy-scan experiments SND and CMD-2, the fit function f (x i ) is the FSR-inclusive bare cross section at the given centerof-mass squared energies s i of the data points, with the fit function being derived from the JHEP02(2019)006 QCD VFF accounting for the kinematic factors from (2.14) and FSR by means of (2.13). In the case of the radiative return experiments, the provided cross-section measurements should be considered as an average value integrated over each energy bin [142,143]. Since the experiments do not provide a bin center weighted by the experimental distribution within the bin, we take as the fit function the theoretical bare cross section including FSR integrated over the energy bins since this prescription should be closest to a reweighting based on the experimental data themselves. The overall effect, equivalent to shifting the bin center according to the theoretical distribution, is small, but becomes relevant in the vicinity of the ρ-ω interference where the VFF is changing rapidly. Finally, we implement the E L bound by adding a penalty term to the χ 2 function that only contributes if the difference between elastic ππ phase and form factor phase is larger than the central value of the bound: where θ is the Heaviside step function. For σ ∆ max i , we use the uncertainty on the bound due to the cross section ratio r. The variation of the elasticity parameter is treated as a systematic uncertainty.
Since the data compilation [118] only considers the contribution to the cross-section ratio r from I = 1 channels, we do not include the isospin-breaking factor G ω (s) in the bound, i.e. we only constrain the phase of the inelasticity factor by identifying but in any case away from the ω resonance the phase of G ω is tiny. In the fit results, we do not include the data points of the E L bound in the counting of the degrees of freedom, otherwise one might encounter a situation where small shifts in the model function change the number of degrees of freedom. This treatment is further justified by the fact that the contribution of χ 2 E L to the total χ 2 is typically very small.

Fit results
In the following, we discuss different fit strategies by comparing the goodness of the fits to single time-like data sets. We also perform simultaneous fits to a single time-like data set and the space-like data sets. These studies allow us to define an optimal fit strategy that we will use in section 4.4 for fits to a combination of time-like (and space-like) data sets.

Fixed ω mass
In a first step, we fix the mass and width of the ω at the PDG values [132]. For simplicity, we use only one free parameter in the conformal polynomial, i.e. N = 2. Therefore, in total JHEP02(2019)006   we have four fit parameters: the two values of the phase shift, the ρ-ω mixing parameter ω , and one parameter c 2 in the conformal polynomial. In table 2, we show the results of the fits to single time-like data sets. Apart from the fit parameters, we show the value for the two-pion HVP contribution to a µ from the energy region √ s ∈ [0.6, 0.9] GeV (including the FSR contribution according to (2.13)). Although the values for a µ are reasonable, the fit quality in general is very poor. The p-values clearly show that these fits are unacceptable.
This conclusion is most severe for the BESIII data set, for which we find a reduced χ 2 of the order of 10. This behavior can be traced back to the statistical covariance matrix, e.g. the exact same difficulties arise for any kind of global fit function. For instance, the Gounaris-Sakurai (GS) [144] fits presented in [74] are performed using the diagonal errors only, while a fit using the full covariance matrix breaks down in the same way as a fit using our dispersive representation. Moreover, this observation stands in marked contrast say to the BaBar data, for which a GS fit was performed including both systematic and off-diagonal statistical uncertainties as well, leading to a χ 2 in a similar range as ours [72]. We conclude that with the statistical covariance matrix as provided together with [74] no statistically meaningful description of the data is possible, and will therefore not consider the BESIII data set in the following. 7

Fitting the ω mass
In a next step, we use the ω mass as a free fit parameter and disregard the input from the PDG. The results of the fits to single e + e − data sets are shown in table 3. The fits to the energy-scan experiments and BaBar are now of good quality. Unfortunately, the fit to KLOE is only improved slightly, and fitting the ω width as well does not improve the fit either.
However, the fit result for the ω mass is not in agreement with the value (4.1) from the PDG [132], which is dominated by e + e − → 3π and e + e − → π 0 γ experiments at SND and JHEP02(2019)006  CMD-2 [64,146,147] as well as frompp → ωπ 0 π 0 [148]. Due to the fact that in the twopion channel the ω resonance appears very close to the broader ρ resonance and only due to an isospin-violating effect, it seems unlikely that the extraction of the ω mass from the two-pion channel should be more reliable than from these channels. Therefore, one might suspect that consistency among the different channels could require allowing for another source of uncertainty related to the energy calibration in the respective experiments. This possible explanation is further pursued below in terms of fits that implement precisely such an energy rescaling. Finally, we note that the fit values of the phase shift are in all fits perfectly compatible with the values (3.1) used in the Roy-equation analysis [10], and, even more importantly, consistent among the different data sets at a level well below the uncertainties quoted in (3.1). This shows the potential in further improving the ππ P -wave phase shift from the present VFF fits.

Energy rescaling
Instead of fitting the ω mass, proper alignment with its PDG value could be ensured by rescaling the energies of the time-like data points i by where ξ j is a small rescaling factor for each experiment j and we have chosen this mapping to leave the two-pion threshold invariant. The rescaling of the energy affects the relation (2.14) between the form factor and the bare cross section (we neglect the rescaling in the FSR correction). The effect can be described by where A(s) ∈ [−1.5, 2] for s ≥ 4M 2 π . The results of this fit are shown in table 4. They are almost indistinguishable from the ones where the ω mass is fit. We also note that the exact form of the rescaling (4.10) proves immaterial, given that a simpler rescaling s i → ξ 2 j s i or a small energy shift are possible as well and lead to almost identical results. As the energy rescaling is at the permille level, the effect on the integrated a µ is very small, while the improvement in the χ 2 compared to the fit with fixed ω mass and no energy rescaling is critical to obtain acceptable fits. In the end, it appears to be simply related to the correct alignment of the ω resonance, which, when insisting on the PDG value (4.1), necessitates some rescaling as in (4.10).  Table 5. The same as table 2, but with both free ω mass and energy rescaling.  However, the implied energy calibration uncertainties as large as 1 MeV in the ρ peak are significantly larger than estimated in the respective experiments [139,142,143], pointing to a corresponding tension in the ω mass among different channels. To separate these issues, in particular in the combined fits, we will from now on keep both a global ω mass and individual rescalings for each experiment as free fit parameters, but constrain each rescaling by an additional penalty ∆χ 2 j = (ξ j M ρ /∆E j ) 2 , where the calibration uncertainty in the ρ peak is estimated as ∆E = 0.4 MeV for the Novosibirsk data sets [139], 0.16 MeV for BaBar [72], and 0.2 MeV for KLOE [142]. In contrast to the E L bound, we will count these terms as additional data points in the number of degrees of freedom. The results shown in table 5 illustrate the fact that a free ω mass and an energy rescaling are all but equivalent, with the corresponding flat direction broken by the requirement that the energy rescaling not be larger than acceptable given the estimate of the experimental calibration uncertainty.

JHEP02(2019)006
In the case of KLOE, we have used the combination of the KLOE08, KLOE10, and KLOE12 results [75], but in all fit variants considered so far the fit quality is significantly worse than for the other experiments. However, as shown in table 6, we observe a significant improvement of the χ 2 if we allow for different energy rescalings ξ j for each of the three KLOE experiments. Since it may indeed be that energy calibration uncertainties differ among the three KLOE data sets, we will allow for three individual rescalings in the following, each constrained by an uncertainty of 0.2 MeV at the ρ peak.

Possible outliers
Let us scrutinize the fit results to KLOE and BESIII. By considering the individual contributions to the χ 2 from each energy bin, we were able to identify in the KLOE set two bins with wildly disproportionate contributions to the χ 2 : if we remove bins #15 and #22 from the KLOE08 set, the total χ 2 reduces by more than 30 units, as shown in table 6 (the set with deleted outliers is marked as KLOE08 ). In figure 6a,  bin contributions, to the χ 2 in the KLOE fit, restricted to the 60 bins of KLOE08. The bins #15 and #22 are marked by red circles. They correspond to 703.56 MeV and 751.66 MeV, where the form factor is expected to show no particular structure. In fact, even in the vicinity of the ρ-ω interference, where the VFF varies much more rapidly, no conspicuous contributions to the χ 2 arise, see figure 6a, while bins #15 and #22 are clearly visible. This suggests to discard them as obvious outliers. We will denote the corresponding data set by KLOE in the remainder of the paper, but also show results for the full KLOE data set. In the end, the main impact is restricted to the goodness of the fit, the results for the HVP contribution to a µ or the pion radius are hardly affected. Unfortunately, in the case of BESIII we were not able to identify similar outliers. In figure 6b, we show the individual bin contributions to the χ 2 in the BESIII fit. We observe fluctuations between huge positive and negative values. If we only take into account the diagonal elements of the covariance matrix, a perfect fit (with a reduced χ 2 around 1) is possible. As mentioned above, this suggests that there might be a problem with the BESIII covariance matrix. We remark in addition that the diagonal elements of the statistical covariance matrix in the supplementary material of [74] do not agree with the diagonal errors published in the same reference.

Including space-like data sets
We now perform fits to a combination of one time-like data set and the space-like data from NA7 (we also tried including F2, in addition, but the gain in statistics is entirely JHEP02(2019)006   Table 7. Combined fits to one time-like experiment and the space-like NA7 data set. In the KLOE set, the two outliers in KLOE08 have been removed. No rescaling of s has been applied to the space-like data.
negligible and we drop the corresponding data set for simplicity). Although for (g − 2) µ we only integrate over the time-like region above the threshold, s ≥ 4M 2 π , analyticity provides the connection between the time-like and space-like region, so that we can use experimental input from both regions to constrain the form factor. In principle, the same discussion of radiative corrections as in section 2.1 arises, but fortunately the applied radiative corrections in the space-like data sets include VP [149,150], so that the provided data for the form factor can be used without further adjustments.
In table 7 we show the results of the combined fits including NA7. The NA7 data are perfectly compatible with the fits to all time-like experiments. Since they have much larger uncertainties than the e + e − experiments, their influence on the fit result is minor, mainly leading to a smaller χ 2 /dof. This is even more so in the case of the F2 data, which do not have any observable influence on the fit.

Varying the order of the conformal polynomial
In a final step, we vary the order N of the conformal polynomial used to describe inelasticity effects. Due to the P -wave constraint (2.23), the number of free parameters in the conformal polynomial is N −1. The fit results for N −1 = 1 . . . 5 are shown in table 8. The fit quality is good in all cases, provided that we remove the two outliers from the KLOE08 set. For small N , the E L bound is fulfilled either automatically or imposed at only one point, while for larger N , the number of points where the bound is activated increases. We have performed fits with up to N −1 = 7 free parameters in the conformal polynomial. In the case of BaBar and KLOE, the χ 2 does not improve any more for N − 1 > 4, while for SND and CMD-2 some further improvement might be inferred, but due to the large number of parameters their fit values become unnaturally large and highly correlated. In all cases, the results for a ππ µ remain stable for larger values of N , with the main effect that the parameters of the conformal polynomial receive large uncertainties and the E L bound becomes an increasingly important constraint. Moreover, the phase of the inelasticity contribution G N in starts to oscillate for higher values of N , indicating further that very large values of N do not correspond to a physically acceptable solution. Therefore, we choose N − 1 = 4 JHEP02(2019)006

Combining data sets
We now present the results of our final fit configuration. We use N − 1 = 4 free parameters in the conformal polynomial, let the ω mass free, and use an energy rescaling for each of the time-like data sets constrained by the expected energy calibration uncertainty, in the case of KLOE three separate rescaling parameters for KLOE08/10/12. From the KLOE08 data set, we remove the two outliers ("KLOE "), but also show the fits to the full set to demonstrate that while the improvement in the χ 2 is significant, the impact on HVP is very minor. All sources of systematic uncertainties described in section 4.1 are considered, leading to the fit results for the parameters M ω , ξ j , ω , and the values of the phase shift at s 0 and s 1 as shown in table 9. The fit errors are inflated by a scale factor S = χ 2 /dof, (4.13) according to the PDG averaging prescription [132]. Next, we perform simultaneous fits to combinations of the data sets. As the fit quality is equally good in all fits to single experiments, we do not introduce any weighting factors, but only apply the inflation factor S (4.13), which increases the fit errors by 12% to 19%. 8 The results of these fits are given in table 10.
In figure 8, we show the fit result for the VFF both in the time-and space-like region together with all the data sets used in the fit. At this scale, the uncertainties of the fit result are barely visible. In figure 9, we show the space-like region of the VFF together with the NA7 data. In figure 10, we focus on the ρ-ω interference region, in order to make JHEP02(2019)006   Table 9. Final fits to single e + e − experiments with N − 1 = 4 free parameters in the conformal polynomial. The first error is the fit uncertainty, inflated by χ 2 /dof, the second error is the combination of all systematic uncertainties. it possible to distinguish between the dense time-like data sets. For comparison, we show in figure 11 the same plot without the energy rescaling (4.10) and for PDG ω mass, so that the effect of the exact alignment of the ω resonance position becomes apparent. In figure 12, we show for the region [0.6, 0.9] GeV the relative deviation of the data points from the fit result, normalized to the fit value of |F V π (s)| 2 . In this variant, one can clearly observe the well-known tension between the BaBar and KLOE data sets [80,81]: the BaBar data lie systematically above the KLOE results, and the fit finds the average as dictated by the experimental covariance matrices.

Extraction of the ω mass
At first sight, it is surprising that our final number for the ω mass resulting from the fit to the combination of all experiments M ω = 781.68 (9)(3) MeV, (4.14) see table 10, is slightly below a naive weighted average of the results from the fits to single experiments in table 9 and even below the fit results to BaBar or KLOE alone. However, most of this effect is explained by the correlations with the other fit parameters. If one performs a multivariate weighted average of the fit parameters where p j is the vector of parameters fit to experiment j and σ j denotes the covariance matrix of the fit parameters p j , one obtains a result very close to the outcome of a fit to the combination of all experiments, in particular, one finds M aver ω = 781.67(9) MeV. A rather large correlation of −39% is present between M ω and ω in the fit to KLOE. This correlation partly explains why the combined value (4.14) lies below the fit results to either BaBar or KLOE, as illustrated in the plot of the error ellipses in figure 13. Small differences between the multivariate weighted average and the result of the combined fit Fit result for the VFF |F V π (s)| 2 Total error Fit error NA7 SND CMD-2 BaBar KLOE08 KLOE10 KLOE12 Fit result for the VFF |F V π (s)| 2 Total error Fit error NA7 VFF fit result and data with energy rescaling Total error Fit error SND CMD-2 BaBar KLOE08 KLOE10 KLOE12 Figure 10. Fit result for the pion VFF in the ρ-ω interference region, together with the e + e − data sets. The data points are shown with the energy rescaling (4.10) and the curve is the fit result with (4.14) for the ω mass.
VFF fit result with M PDG ω and data without energy rescaling Total error Fit error SND CMD-2 BaBar KLOE08 KLOE10 KLOE12 Figure 11. Fit result for the pion VFF in the ρ-ω interference region, together with the e + e − data sets. The curve is the result of the VFF fit to the data points including energy rescaling as shown in figure 10, but with an ω mass reset to the PDG value and compared to the original data points without energy rescaling.
Relative difference between data sets and fit result Total error Fit error SND CMD-2 BaBar KLOE08 KLOE10 KLOE12 Figure 12. Relative difference between the data points (including the energy rescaling (4.10)) and the fit result for the VFF, normalized to the fit result for |F V π (s)| 2 . As in all plots, we show fit errors and total uncertainties as two separate error bands. The total uncertainty is given by the fit error and the systematic uncertainty, added in quadrature. All e + e − (KLOE ), NA7 SND CMD-2 BaBar KLOE Figure 13. Error ellipses for the parameters ω and M ω resulting from fits to single experiments and the fit to the combination of all experiments. The smaller ellipses are standard error ellipses that correspond to ∆χ 2 = 1, the larger ellipses are inflated by the scale factor (4.13).   Table 10. Final fits to combinations of data sets with N − 1 = 4 free parameters in the conformal polynomial. The first error is the fit uncertainty, inflated by χ 2 /dof, the second error is the combination of all systematic uncertainties, apart from the case of the parameters c i in the second panel, where it includes only the uncertainties due to the phase input and Γ ω , but not the variation in s c and N . The third panel gives the energy rescalings that belong to the respective fits.

JHEP02(2019)006
can be observed for the other parameters, which is a sign of the non-linear dependence of the fit function on the parameters.
Both our final result (4.14) for the ω mass and the results from fits to single experiments in table 9 disagree with the PDG average (4.1). Since this discrepancy is not driven by a single experiment, it seems to indicate that e + e − → π + π − data indeed unanimously favor an ω mass substantially lower than the current PDG average. The latter is mainly based on M ω = 782.79(8)(9) [146], M ω = 782.68(9)(4) [64] from e + e − → 3π, but also includes M ω = 783.20(13)(16) [147] from e + e − → π 0 γ and M ω = 781.96(13)(17) MeV from pp → ωπ 0 π 0 [148], both of which do not influence the average much (or at least cancel each other effectively). Our determination is thus closer to thepp reaction, but in conflict with the 3π data.
This observation is not new, see e.g. [39], and the mismatch has already been pointed out by the BaBar collaboration [72]. The latter analysis is worthwhile reviewing in some detail here because it may offer an indication which effects could be contributing to the puzzle. In [72], the data are analyzed with the help of a sum of GS parametrizations for all the relevant resonances in their energy range, including relative phases between the various GS terms, which are allowed to float and determined by a fit to the data. In the unconstrained fit they obtain M ω = 781.91 (18) (16) MeV and a relative phase φ ρω = −0.6(2.1) • ,

JHEP02(2019)006
compatible with zero. This is in very good agreement with (4.14), but not with the PDG. They then observe that there is a strong correlation between φ ρω and M ω and that if one fixes the ω mass at the PDG value one can still obtain a good fit: in this case, however, the value of φ ρω changes to 7.8(1.3) • , which is compatible with the phase obtained by CMD-2 in their fit to the e + e − → π + π − data [68]: φ ρω = 10.4(3.8) • . Conversely, if they constrain this phase in the fit to the BaBar data to the CMD-2 value, the ω mass becomes M ω = 782.68 (12) (27) MeV, perfectly compatible with the PDG.
Assuming that the fit quality has not decreased significantly in the latter fit, we can conclude that if CMD-2 and BaBar use similar parametrizations to fit their respective data, the values of the parameters they obtain are compatible with each other (taking into account correlations). Still, the first unconstrained fit shows that the BaBar data prefer a value of the ρ-ω phase compatible with zero, which is a very good sign because the presence of such a phase would violate analyticity and unitarity (as it would give a complex form factor even below the ππ threshold). In our framework such a phase is strictly forbidden and the agreement with the unconstrained BaBar fit both on the absence of this phase as well as on the value of the omega mass is very reassuring.
This raises the question whether also other determinations of the omega mass have used unphysical parametrizations, and the answer is unfortunately positive: the determination based on e + e − → π 0 γ [147] includes such a relative phase of 13.3 • , which might explain the resulting value for M ω even higher than from the 3π channel. In contrast, the extractions based on e + e − → 3π do not include a Breit-Wigner ρ-resonance in their parametrization, only the ω and φ resonances together with a smooth background, but the complex phases between ω, φ, and background could still distort the extracted ω mass. Both for e + e − → 3π and e + e − → π 0 γ representations exist that do not suffer from these shortcomings [53,56]. In these papers, good fits were obtained while using the PDG ω parameters as input, which indicates a substantial model dependence in the extraction from e + e − → π 0 γ [147], but likely only a small effect in e + e − → 3π [64,146]. For a firm conclusion more thorough fits to 3π and π 0 γ data including the respective uncertainties would be necessary.
For these reasons, the high significance of the discrepancy, more than 5σ if taken at face value, is puzzling, in particular given that the extraction from the isospin-conserving 3π channel should, in principle, be more reliable than the isospin-breaking effect in e + e − → π + π − . Another potential subtlety could concern the definition of the ω mass in view of electromagnetic corrections, but estimates of the corresponding effect [114] ∆M ω = e 2 2g 2  Table 11. Values for a ππ µ from our final fits to single e + e − experiments. The first error is the fit uncertainty, inflated by χ 2 /dof, the second error the combination of all systematic uncertainties. We provide the results for several energy regions separately, to enable a detailed comparison with other (future) evaluations. The energy regions in the third block are provided to facilitate comparison with [37] and the results of the direct integration [75].
In this paper, we aim to derive the HVP contribution to a µ based on the available experimental information on e + e − → π + π − subject to a comprehensive analysis of the constraints from analyticity and unitarity. From this point of view, there is no indication to assume a common systematic effect in all experiments that would restore agreement with the 3π channel, we will therefore pursue the analysis of the HVP contribution based on the fits in the preceding subsection. However, in addition to the known tension between the KLOE and BaBar data, this discrepancy in the ω mass extracted from the 2π and 3π channels deserves further attention.

Consequences for the anomalous magnetic moment of the muon
In  Table 12. Values for a ππ µ from our final fits to combinations of data sets. The first error is the fit uncertainty, inflated by χ 2 /dof, the second error is the combination of all systematic uncertainties. the result for the combined fit does not exactly coincide with a naive weighted average of the fit results to single experiments: most importantly, correlations play a role in the same way as discussed in section 4.5. Further small deviations are due to the non-linear dependence of the fit function on the parameters, which leads to distortions of the χ 2 . We checked that the deviations of the χ 2 from a quadratic function in the parameters are very small within the standard confidence regions of the parameter space. Further away from the χ 2 minimum, these deviations become more important and they have an observable effect in the combination of the BaBar and KLOE data sets, which reflects the well-known tension between these two experiments, see figure 12. Taking into account the correlation of the systematic uncertainties, this discrepancy between the BaBar and KLOE results for a ππ µ below 1 GeV amounts to 2.6σ. Finally, the relative size of various sources of systematic uncertainties is illustrated in figure 15. The dominant systematic error is due to the order of the conformal polynomial, followed by the Roy parameters (including ι 1 ) and s c from the conformal expansion. This pattern holds true for most of the fit variants considered.
Where published results are available, we have included the comparison in the tables, e.g. from direct integration [75,81] and the dispersive analysis [37]. We find that in those cases where reference values exist, our results appear well compatible, within uncertainties of a similar size. An exception is the comparison to the direct integration of the data JHEP02(2019)006 between √ 0.1 and √ 0.95 GeV performed by KLOE [75] where our method shows a significant reduction of the uncertainties: this is mainly due to the region below 0.6 GeV where KLOE data show a loss of precision. With our approach once precise data are available in the most sensitive region around the ρ peak they strongly constrain the curve in the whole low-energy region and the extrapolation down to the two-pion threshold does not lead to any loss in precision: this is a clear advantage of our method with respect to the application of the trapezoidal rule. On the other hand, tables 11 and 12 show that in the regions where there are high-quality data, these are so precise and densely spaced that our method does not lead to an increase of precision, but mainly serves as a check of the consistency of the data with the principles of analyticity and unitarity. Reversing the argument, the fact that our uncertainties are of similar size as those of other analyses shows that the systematic uncertainties in the dispersive representation, which we have investigated in detail, are well under control in the whole region below 1 GeV. We stress that our uncertainty estimates, illustrated and summarized in figure 15, rely on minimal assumptions, the dispersive parametrization as a consequence of QCD and the covariances matrices provided by experiment, where the latter then effectively determine the relative weight of each data set in the combined fit. In particular, a local inflation of the uncertainties would be difficult to justify in this formalism, which emphasizes the importance of the finding that each data set allows for a statistically acceptable fit once potential uncertainties in the energy calibration are taken into account (and the two outliers in KLOE08 removed). 9 We look forward to more detailed comparisons to direct integration [80,81], which should lead to a better understanding of the uncertainties in the critical ππ channel and thereby to a consolidation of the overall uncertainty estimate for HVP.
Our most comprehensive result gives the full contribution below 1 GeV in a combination of all available time-and space-like constraints where the inclusion of the space-like data does allow for a modest reduction of the uncertainty from 2.8 to 2.6 units. As noted before [37], the main advantage over direct integration occurs in energy regions where data are still scarce, most notably the low-energy region Our result agrees with the combination of e + e − data sets from [37], a ππ µ | ≤0.63 GeV = 133.0(8) × 10 −10 , which provides another important cross check given several conceptual differences compared to our study. 10 The main difference to our approach concerns the fact that the ππ phase shift is not fit to the data, but taken as an input. The dispersive 9 As demonstrated by table 12, the central values in the combined fit to all experiments barely change when the two KLOE08 outliers are retained: the value (5.1) for a ππ µ below 1 GeV increases by 0.2 × 10 −10 , but the total χ 2 is worse by about 30 units and leads to a slightly larger scale factor (4.13) of 1.13 instead of 1.11.
10 Note that the final number quoted in [37], a ππ µ | ≤0.63 GeV = 133.3(7) × 10 −10 , also includes information from τ data, but given the difficulties in controlling the required isospin-breaking corrections we only consider e + e − data here. formalism is then set up in such a way that the phase shift in the elastic region alone, in combination with data for the modulus of the VFF in the energy region √ s ∈ [0.65, 0.71] GeV, constrains the VFF in the low-energy region √ s ≤ 0.63 GeV. On the one hand, in this way the systematic uncertainties related to the high-energy continuation of the phase shift and the inelastic corrections no longer need to be considered, but on the other hand the method is then necessarily restricted to rather low energies. In contrast, our representation remains applicable as long as inelastic corrections can still be controlled, within the formalism that we have employed here at least up to 1 GeV. Moreover, our approach avoids a circularity problem that arises because the ππ phase shifts used as input have been extracted from previous form factor fits themselves, even though the numerical impact of this effect might be negligible in the end. The HVP result for the low-energy region agrees in both implementations of dispersive constraints on the pion VFF. Both phase shift values are fully compatible with the ranges from (3.1), with statistical uncertainties well below these errors. In all cases, we observe that the fit results are extremely stable among different data sets, in such a way that by far the dominant uncertainty now arises from systematic effects.
To arrive at (6.1), we considered separately each of the 25 additional parameters in the Roy solution, see section 3, and added all uncertainties in quadrature only at the very end of the calculation. However, very similar results emerge if instead one defines a smooth band around the central Roy solution by adding in quadrature all uncertainties other than those from δ 1 1 (s 0 ) and δ 1 1 (s 1 ). The propagation of the individual parameter uncertainties also allows one to identify the source of the relatively large systematic effects in δ 1 1 (s 1 ), which are dominated by the asymptotics of the imaginary part of the partial wave, Im t 1 1 , as well as a low-energy parameter from the isospin-0 S-wave. This interrelation shows that for a global analysis of low-energy ππ phase shifts the role of these systematic effects, in particular the interplay with the Roy parameters corresponding to other isospin channels, needs to be carefully investigated. This will be addressed in future work.
In this regard, it might appear curious that the final error quoted for δ 1 1 (s 1 ) is actually slightly larger than in (3.1). However, one should keep in mind that in the solution of the Roy equations [10], all the parameters are to be varied independently within their uncertainty ranges. With our fit of the VFF to data, the phase values (6.1) are no longer independent parameters but correlated with the remaining Roy parameters p i . Linearizing the fit result around their central values p c i , we write in order to make the systematic dependence on the additional Roy solution parameters explicit. The values ofδ 1 1 (s i ) now only contain the systematic effects that are independent of the 25 additional parameters of the Roy solution: Fit result for the ππ P -wave phase shift δ 1 1 Systematic error Fit error Figure 16. Fit result for the elastic P -wave ππ scattering phase shift δ 1 1 . The gray band shows the systematic uncertainty due to the parameters in the dispersive form factor representation, while the black error band representing the fit uncertainty is hardly visible. Note that, in contrast to figure 3, the band includes the systematic uncertainties related to the asymptotic continuation of the phase shift. 25 remaining Roy parameters. In particular, this separation clearly shows the improvement in the determination of the phase shift compared to the independent parameter ranges (3.1).
As illustrated by figure 15, these issues are immaterial for the HVP contribution, so that for the present application we do not attempt to reduce the systematic errors further. The present status of the P -wave phase shift, corresponding to (6.1), is illustrated in figure 16. As expected, the band characterizing the systematic uncertainties widens rapidly above 1.15 GeV, while throughout the contribution of the statistical error is completely negligible.

Charge radius of the pion
The charge radius of the pion, r 2 π , is defined by the derivative of the VFF at s = 0 where the derivative is again evaluated via a dispersion relation. With the VFF determined from the fit to e + e − → π + π − data, this integral produces the results collected in table 13.  Table 13. Charge radius corresponding to the fits to single time-like experiments and to combinations of data sets. The errors in the first five columns are the fit uncertainties, inflated by χ 2 /dof. The results in the last column correspond to N − 1 = 1. The first error is the inflated fit uncertainty, the second error is the total uncertainty (which includes the variation N − 1 = 1 . . . 5). extracted value for r 2 π in N , to the extent that the most stable results are obtained for small values of N and, as seen from table 13, N − 1 = 4 and 5 already begin to go astray. We still keep the full systematic variations for N − 1 = 1 . . . 5, otherwise one would have to investigate in more detail the potential role of inelastic effects above the energy range constrained by the E L bound. As central values we quote the results for N − 1 = 1, both motivated by the fact that the extrapolation uncertainties of the conformal polynomial beyond 1.15 GeV are expected to be smallest for the lowest order and since these values happen to lie around the middle of the range given in table 13. 11 Our final result, including both time-and space-like data sets, reads r 2 π = 0.429(1)(4) fm 2 = 0.429(4) fm 2 . (7.2) Within uncertainties, this value is consistent with the previous dispersive extraction r 2 π = 0.432(4) fm 2 from [151], but the tension with the PDG average r 2 π = 0.452(11) fm 2 [132] is further exacerbated. However, as noted before [39], 12 this average does not contain any modern e + e − → π + π − data sets and, if potentially model-dependent extractions from eN → eπN [152,153] were excluded, would be dominated by NA7 r 2 π = 0.439(8) fm 2 [134], in better agreement with (7.2). Indeed, if the NA7 data were in conflict with our dispersive determination, a simultaneous fit of time-and space-like data would not be possible. Our calculation therefore provides further evidence that the PDG average for r 2 π needs to be revised. 11 For a central value defined by N − 1 = 4, the final result would change to r 2 π = 0.426(1)(6) fm 2 , where the systematic error points entirely in the upward direction. 12 We observe that the results in table 13 do not change within the fit uncertainty of 0.001 fm 2 if VP is absorbed into the definition of the VFF, whereas significant effects do occur in the evaluation of a ππ µ .

JHEP02(2019)006 8 Summary and outlook
Analyticity and unitarity imply strong constraints both on ππ scattering and the pion VFF.
In this paper, we analyzed these constraints comprehensively as regards consequences for the HVP contribution to the anomalous magnetic moment of the muon, including both a consistent implementation of the experimental uncertainties as well as the systematic uncertainties associated with the dispersive representation. The central outcome of this study (5.1) a ππ µ | ≤1 GeV = 495.0(1.5)(2.1) × 10 −10 = 495.0(2.6) × 10 −10 , shows that the main complications in such a representation, arising from inelastic corrections and high-energy contributions, can be controlled at a level that renders the dispersive approach a valuable complementary perspective to the direct integration of the experimental data. In particular, it provides the best controlled extrapolation down to the two-pion threshold where data are less precise or just absent.
With the present analysis we have therefore laid the ground work to consolidate the uncertainty estimate for the ππ channel in HVP. In contrast to the direct integration, we cannot allow for the local inflation of uncertainties since the dispersive fit function defines a global constraint. For that reason it is critical that once possible uncertainties in the energy calibration are taken into account all present data sets can be described in a statistically acceptable way, providing a strong check on their internal consistency. The combination of data sets then follows in a straightforward way from the propagation of the uncertainties incorporated in the covariance matrices provided by experiment, up to a small inflation of the final uncertainties by χ 2 /dof ∼ 1.1 in the standard manner, a global scale factor that is much smaller than the local scale factors up to 3 that are required otherwise. In this way, we have obtained a combination of the available e + e − data sets with minimal assumptions, relying only on the global fit function that follows from QCD and the stated experimental uncertainties. We expect that a future more detailed comparison with direct integration should lead to a better understanding of the uncertainties in the ππ channel and eventually make the overall error estimate more robust. As an added benefit, a dispersive approach has to be able to accommodate space-like data sets at the same time, which not only provides a further consistency check both on the data and the formalism, but in this case actually leads to a modest reduction in uncertainty.
At this point, the experimental data on e + e − → π + π − are so precise that the systematics of the dispersive representation begin to dominate, an observation that becomes most apparent for the values of the ππ phase shift extracted from the fit (6.1) and (6.4). For HVP the mismatch between statistical and systematic errors is still relatively small, but for future data sets improved variants of the dispersive representation could be investigated. For instance, figure 15 shows that by far the dominant effect arises from the order N of the conformal polynomial that describes the inelastic corrections, which we estimated very conservatively by the maximum deviation found among all statistically meaningful fits.
Here, more precise data, in combination with the E L bound, might allow one to actually identify an optimal value or range of N or even attempt an explicit description within the dispersive approach of inelastic effects in terms of physical processes and thereby signifi-

JHEP02(2019)006
cantly reduce the associated uncertainty. Another issue that warrants further investigation concerns the mass of the ω, for which it would be important to clarify the origin of the current mismatch between extractions from the 2π and 3π channels.
Beyond the HVP contribution, our results for the VFF are important for an improved understanding of low-energy ππ scattering. While most recent work has focused on improving the isospin-0 S-wave, the input used for the isospin-1 P -wave in solving the full system intertwined by crossing symmetry actually goes back to by now outdated analyses of the pion VFF. Based upon the present work it will become possible to perform a global analysis of low-energy ππ phase shifts including the stringent constraints on the P -wave from the modern high-statistics e + e − → π + π − experiments.