Test of analyticity and unitarity for the pion form-factor data around the $\rho$ resonance

High-statistics data on the $e^+e^-\to \pi^+\pi^-$ cross section and the pion vector form factor have been obtained recently by several collaborations. Unfortunately, there are some tensions between different datasets, especially the most precise ones, which have not been resolved so far. Additional independent constraints on the data are therefore of interest. We consider a parametrization-free method of analytic extrapolation proposed recently, which is based on a mixed phase and modulus extremal problem and combines rigorous upper and lower bounds with numerical simulations to account for the statistical distributions of the input and output values. Spacelike data on the form factor and measurements of the modulus in the region $(0.65-0.71)$ GeV are used as input. In previous works, the formalism was applied for extrapolating the form factor to low energies. In the present work, we use it as a stringent and model-independent test of consistency with analyticity and unitarity for the high-statistics data around the $\rho$ resonance. The study reveals some inconsistencies, in particular below the $\rho$ peak the BABAR data are slightly higher than the band of extrapolated values, while above the $\rho$ peak all the data are situated at the lower edge of the band. The implications of the results on the two-pion vacuuum polarization contribution to the anomalous magnetic moment of the muon are briefly discussed.


I. INTRODUCTION
The electromagnetic form factor F V π (t) of the pion is a fundamental quantity for strong-interaction physics, intensively studied theoretically and experimentally for more than sixty years. The recent interest in its precise determination is mostly driven by the anomalous magnetic moment of the muon, a µ = (g − 2) µ /2. There is at present a disagreement between the experimental value of the muon anomaly measured at BNL [1] and its evaluation in the Standard Model (SM). Two new experiments aim at reducing the experimental uncertainty by a factor of four: the E989 experiment at Fermilab [2], which started running in 2018, and the E34 experiment at J-PARC, which plans to start its first run in 2024 [3]. In parallel, there is a continuous effort for improving the accuracy of the theoretical calculation of a µ in the SM (for a recent review and earlier references see [4]).
The hadronic vacuum polarization (HVP) contribution to the muon g − 2 is the leading hadronic contribution, which cannot be calculated using perturbative QCD and dominates the theoretical uncertainty. By analyticity and unitarity, the hadronic loops contribute to leading order (LO) as a dispersion integral over the cross section of the e + e − annihilation into hadrons. This allows a datadriven evaluation of HVP using experimental data (the most recent determinations are reported in [4][5][6]). At low energies, below √ s = 1 GeV, the dispersive integral is dominated by the cross section of the e + e − annihilation into two pions. Several high-statistics e + e − experiments [7]- [15] have been designed and operated recently in order to measure this quantity with increased precision. However, as noted in [4][5][6], there are some tensions between the data of different experiments, particularly the most precise ones, BABAR [10,11] and KLOE [12][13][14].
It is useful to keep in mind that to LO the cross section of the e + e − → π + π − process is expressed in terms of the modulus squared of the pion vector form factor, F V π (t). This allows, in principle, to use additional theoretical and experimental information available on the form factor in order to improve the accuracy, especially in regions where the data are still poor. The most powerful approach for the study of the form factor is dispersion theory, which exploits analyticity, unitarity and crossing symmetry. The form factor has been calculated also on the lattice, but the results did not reach the same accuracy until now, and we shall not consider them.
The dispersion theory of the pion form factor traditionally exploits Fermi-Watson theorem [16,17], which states that below the first inelastic threshold the phase of the form factor on the unitarity cut is equal to the P -wave phase shift of ππ elastic scattering. By the wellknown Omnès representation, the form factor is reconstructed as an analytic function in the whole complex t-plane from its phase of the boundary. The most recent dispersive analysis [18] is based on a parametrization involving an Omnès function multiplied by factors which account for inelastic and isospin-breaking effects. The P -wave ππ phase shift used as input was obtained by solving Roy equations, which fully exploit analyticity, unitarity and crossing symmetry of pion-pion scattering amplitudes, and the free parameters of the model have been fixed by fitting data on the modulus of the form factor from e + e − experiments.
As noted in [18], in the combined fit of all data one can observe the well-known tension between the BABAR and KLOE data: the BABAR data lie systematically above the KLOE data. The fit follows the average of the two in accordance with their respective covariance matrices. In the Omnès representation used in [18], an assumption about the phase in the inelastic region, where it is not known, has been adopted. Actually, the parametrization contains an additional factor, which is complex above the inelastic threshold and can modify the initial phase. However, the independence of the results on the choice of the phase in the Omnès function is not formally proved.
As remarked for the first time in [19], in order to avoid assumptions on the unknown phase in the inelastic region, one can use instead the phase some information available in this region on the modulus. This leads to a mixed phase and modulus representation, which can be generalized by including also an arbitrary number of values of the form factor at points inside the analyticity domain. The price to be paid for the independence on the unknown phase above the inelastic threshold is that the formalism cannot predict definite values, but only upper and lower bounds on the form factor. The bounds can be calculated exactly with techniques of functional analysis in terms of the information used as input (for technical detail see [20,21]).
An important step forward was achieved in Refs. [22][23][24][25], where the rigorous bounds derived from analyticity have been merged with numerical simulations which properly take into account the statistical distribution of the input and output values. This elaborate formalism has been applied for the determination of the pion charge radius [22] and the modulus of the form factor at low energies [23][24][25]. In particular, a more accurate value of the HVP contribution to a µ from energies below 0.63 GeV has been obtained. The input used in these works consisted from the most recent measurements of the form factor in the spacelike region and the data on the modulus measured in the range of (0.65 − 0.71) GeV, where the various experiments are in somewhat better mutual agreement than in other regions.
In principle, the method can be used to predict the form factor also at higher energies, above 0.71 GeV, up to the first inelastic threshold. However, as will be clear below, the analyticity bounds become gradually weaker at energies far from the input range, especially when the inelastic threshold is approached. Therefore, unlike the low-energy region where the predictions have been more precise than the experimental data, above the ρresonance region, where precise data by high statistics experiments are available, we do not expect our determinations to compete with them in precision. However, the formalism is still useful since it provides a stringent test of consistency of the form-factor data with analyticity and unitarity.
The purpose of the present work is to explore the outcome of this test for energies below and above the ρ resonance. The paper is organized as follows: in the next section, we give a brief description of the formalism, expressing it as a parametrization-free test of consistency for the values of the form factor at various energies. In Sec. III we describe the theoretical and phenomenolog-ical information used as input and the methodology for extrapolating the form factor in the output region. In Sec. IV we present our results, confronting the extrapolated values with the experimental data in the region (0.72 − 0.9) GeV as the modulus input are taken from the region (0.65 − 0.71) GeV. A discussion of the results, in particular of their implication on the HVP contribution to the muon g − 2, is given in Sec. V. The paper has an Appendix, which contains an explicit proof of the independence of the results on the unknown phase of the form factor above the inelastic threshold.

II. PARAMETRIZATION-FREE ANALYTICITY AND UNITARITY CONSTRAINTS
The pion vector form factor F V π (t) is an analytic function in the complex t plane cut along the real range t ≥ t + , where t + = 4m 2 π . It satisfies the Schwarz reflection property F V π (t * ) = (F V π (t)) * and the normalization F V π (0) = 1. Along the cut the form factor is a complex function. According to Fermi-Watson theorem [16,17], below the first inelastic threshold the phase of the form factor is equal to the P -wave phase shift δ 1 1 (t) of the ππ elastic scattering. Since this theorem is valid in the exact isospin limit, we must first remove the main isospin-violating effect, known to arise from ω − ρ and φ − ρ interference.
We shall follow the standard approach [18,26,27] to do this, by defining a purely I = 1 function F (t) as where the function F ω+φ (t), specified in the next section, includes the I = 0 contribution due to ω and φ resonances. Then Fermi-Watson theorem allows us to write where one can assume with a good approximation that the first inelastic threshold is set up by the ωπ-production threshold, and take √ t in = m ω + m π = 0.917 GeV. In the elastic region, the phase shift δ 1 1 (t) is known with high precision from dispersion theory for ππ scattering [28][29][30].
Above t in , where the phase of the form factor is not known, we use the information available on the modulus from experimental measurements and perturbative QCD. Specifically, we adopt a conservative condition written as where the integral converges and allows an accurate evaluation of I from the available information.
From the phase and modulus conditions (2) and (3), using techniques of functional analysis and optimization theory (for a review see [21]), one can derive constraints on the values of the form factor and its derivatives at points inside the analyticity domain. Omitting the proof given in [20,21] (see also the Appendix of [24]), we shall write down these constraints in a compact form, by expressing the form factor in the isospin limit as where and the functionz conformally maps the t complex plane cut for t ≥ t in onto the unit disk |z| < 1, such thatz(0) = 0 and the upper(lower) edges of the cut become the upper(lower) semicircles in the z plane. We recall that O(t) is an Omnès function, where φ(t) is equal to δ 1 1 (t) for t ≤ t in and is an arbitrary function above t in . The function denoted as ω(t) is analytic in the t plane cut for t > t in and has the modulus equal to |O(t)| on the cut, and w(z) is an outer function (analytic and without zeros in the unit disk |z| < 1), with modulus on the boundary |z| = 1 of the disk related to to the weight in the integral (3).
Finally, the function g(z) appearing in (4) is an analytic function in the disk |z| < 1 and satisfies the boundary condition This condition implies rigorous correlations among the values of the function g(z) and its derivatives at points inside the holomorphy domain, |z| < 1. In particular, if z n ∈ (−1, 1), n = 1, 2, 3 are three arbitrary real points, the following positivity condition holds where D is the determinant defined as where and Moreover, the principal minors of the determinant D must also be nonnegative. We emphasize that the normalization condition F V π (0) = 1 is included as input and is implemented in the definition of g(0) in (12).
The inequality (10) defines an allowed domain for the real values g(z n ), which can be expressed in a straightforward way in terms of the values of the form factor at the points t n =t(z n ), wherẽ is the inverse of (8).
In particular, for t n < t + , when both F (t n ) and O(t n ) are real, we have from (4): while for t + < t n < t in , when F (t n ) and O(t n ) are complex and have equal phases, where the modulus |O(t)| of the Omnès function is obtained from (5) by the principal value (PV) Cauchy integral .
Using these relations, the inequality (11) can be written as an explicit consistency test imposed by anayticity on the values F (t n ) at three arbitrary points below t in . As proved in [20,21], the constraint exploits in an optimal way the input information. Moreover, although the function O(t) depends on the phase φ(t) for t > t in , the test is independent on it. For completeness, we give in the Appendix, along with the formal argument put forward in [19,20], an explicit proof of this independence, presented here for the first time.
In our analysis, we have treated two of the values of F (t) in (11) as input, using experimental information on the form factor (or its modulus) at two points on the spacelike and timelike axes. For each input, the inequality (10) becomes a quadratic constraint on the form factor at the third point, from which upper and lower bounds are exactly derived. The test will consist from the comparison of these bounds with the experimental data available in the third region. The regions used as input and output will be specified below.

A. Inputs
We have extracted the pion form factor from the cross section of the e + e − → π + π − process using standard correction factors (see for instance Appendix B of [23]). We then converted the input values of F V π (t) into the isospinconserving function F (t) using (1), solved the optimization problem for this function and finally re-expressed the results in terms of |F V π (t)|. The function F ω+φ (t), which accounts for the isospin violation due to ω − ρ and φ − ρ mixing, has been taken of the form [26,27]: where = −2 × 10 −3 , = 5 × 10 −3 and the masses and the decay widths are taken as m ω = 781.86, m φ = 1019.46 MeV, Γ ω = 8.5 MeV and Γ φ = 4.3 MeV [31].
We have taken the phase shift δ 1 1 (t) entering the integral (5) below t in from [28,29] (Bern phase) and [30] (Madrid phase). They have been calculated with good precision using dispersion relations for ππ scattering and are mutually consistent, with slightly larger uncertainties of the Bern phase near t in . Above t in we have taken for φ(t) an arbitrary expression. As we have already mentioned, the results are independent of this arbitrariness. The proof of this is provided in Appendix A.
We have calculated the integral (3) using the BABAR data [10] from t in up to √ t = 3 GeV, smoothly continued with a constant value for the modulus in the range 3 GeV ≤ √ t ≤ 20 GeV, and a decrease ∼ 1/t at higher energies, as predicted by perturbative QCD [32,33]. With this input we have obtained I = 0.578 ± 0.022, where the uncertainty is due to the BABAR experimental errors. In the calculations we have used as input for I the central value increased by the error, which leads to the most conservative bounds due to a monotonicity property discussed in [23].
The input at interior points was taken from the most recent experimental measurements of F π Collaboration at JLab [34,35] on the spacelike axis,  The choice of the input region (0.65 − 0.71) GeV was motivated in our previous works [22][23][24] by the fact that here the data have a better precision than at the lower energies where we extrapolated the form factor 1 . The number of input points in this region from each experiment is given in Table 1 of Ref. [23]. We adopt the same input region in the present work, since the determinations of different experiments, especially the most precise ones, BABAR and KLOE, are more consistent among them in this region than at higher energies, around the ρ resonance (see for instance the comparisons in Fig. 13 of [4]).
We recall that the hadronic decays of τ leptons have been used in the past as an alternative source of data in the evaluation of the HVP contribution to a µ . Isospinbreaking (IB) and electromagnetic corrections must be applied in order to convert these data to the cross section of e + e − annihilation (for a recent evaluation and earlier references see [36]). However, as remarked in [4], the present understanding of these corrections is not yet at a level allowing the use of τ data in the muon g − 2 determinations. In view of this consensus in the community, in this work we will use only e + e − data.

B. Methodology of |F V π | determination
We have taken the range (0.72−0.9) GeV as the output region, where we extrapolate the form factor and confront the results with the experimental data.
A nontrivial complication is the fact that the experimental values used as input are known up to statistical and systematic uncertainties. This requires to properly merge the formalism of analytic bounds with statistical simulations. The problem was solved in Refs. [22][23][24] by generating a large sample of pseudodata, achieved by randomly sampling each of the input quantities with specific distributions based on the measured central values and the quoted errors. For each point from an input statistical sample, upper and lower bounds on |F V π (t)| at points t in the output region have been calculated using the formalism described in the previous section. Note that the input points had to pass a consistency condition in order to be included in the sample. Indeed, some of the minors of the determinant D, mentioned below (13), involve only input quantities and, if the positivity condition is violated, the corresponding points had to be rejected. Finally, a set of allowed values in between the bounds has been uniformly generated, taking into account the fact that all the values between the extreme points are equally valid. The number of generated points was adapted to the width of the allowed range, being larger fow wider ranges.
In this way, for a specified spacelike and timelike input, we generated a large sample of output values of |F V π (t)| at each point √ t in range (0.72 − 0.9) GeV of interest. The output distributions turn out to be close to a Gaussian, allowing the extraction of the mean value and the standard deviation (defined as the 68.3% confidence limit interval). The values obtained with input from each experiment at different energies and of the different experiments (SND, CMD2, BABAR, KLOE11, KLOE13 and BESIII) have been then combined using an averaging prescription proposed in [37], where the correlations between different values are derived from the values themselves. The two spacelike input points have been considered separately and then we took the average of the two predictions for both the central value and the error. The same conservative procedure has been adopted for the phases: we considered the results obtained separately with the Bern and Madrid phases and also the simple average of the corresponding predictions. The comparison of these results with the experimental values serves as a test of unitarity and analyticity of the experimental determinations in the input and output regions.

IV. RESULTS
We start by illustrating the statistical distributions of the output values of |F V π (t)|, obtained with a specific input. In Figs. 1 we show these distributions at several values of t, obtained with the pseudodata sample generated with the experimental measurement at 0.699 GeV by BABAR, the first spacelike input (19) and Bern phase. One can note the increasing number of points in the samples with increasing t. The explanation is that the allowed bands between the upper and lower bounds become wider at larger energies, and more intermediate points must be generated in order to obtain the correct statistical distribution. All the distributions are very close to a Gaussian, allowing the extraction of a central value and a standard deviation at 68.3% confidence level.
From these distributions, by applying the averaging procedure described in Sec. III B, we obtained a central value and a standard deviation for the modulus for all energies in the output region (0.72 − 0.9) GeV. For convenience, we present the extrapolated values of |F V π (t)| 2 at 68.3% CL in this region as a band denoted as "allowed band". This is the main prediction of our method, which we compare with the experimental values available in the same region.
In Fig. 2 we show the allowed bands for |F V π (t)| 2 in the range (0.72 − 0.9) GeV, obtained by using as input all e + e − experiments (SND, CMD2, BABAR, KLOE 11, KLOE 13 and BESIII) in comparison with the e + e − experimental data in the same region. The bands obtained with the Bern and Madrid phases are shown separately.
Several remarks can be made from this figure. First, below the ρ peak the bands of extrapolated values are rather narrow, competing with the experimental data in this region. Only BABAR and KLOE data are more precise than the prediction. Several BABAR points lie definitely above the allowed band, while several KLOE points lie below the band. We can say that these points are in certain tension with analyticity and the data from the lower-energy region, (0.65 − 0.71) GeV. For the other datasets, which have larger uncertainties, no inconsistencies are seen.
Above the ρ peak, the allowed bands obtained by extrapolation become gradually wider, as anticipated, and cannot compete in precision with the experimental data. However, some interesting features can be seen also in this region. All the data lie at the lower edge of the wider band obtained with the Bern phase [29], and are definitely below the band obtained with the Madrid phase [30], which is more narrow since the quoted uncertainties are smaller. Some experimental points near 0.9 GeV are also below the average band obtained with the two phases, presented in Fig. 3. The above predictions have been obtained by using as input all the data in the range (0.65 − 0.71) GeV. It is of interest to see what happens if we restrict the input to some datasets. In Fig. 4 we show the allowed band obtained using as input only BABAR 2009 data [10] and Bern phase. One can see that now the BABAR points below the ρ peak are in more agreement with the allowed band, while the KLOE data lie definitely below the band and above the ρ peak all the experimental data are well below the allowed band. As an opposite choice, we show in Fig. 5 the allowed band obtained using as input all the data except BABAR, with Bern phase used as before. The data are now in more agreement with the allowed band, except the BABAR data below the ρ peak which are definitely higher. Finally, in Figs. 6 and 7 we present the allowed bands obtained using as input only KLOE 2011 [13] and KLOE 2013 [14] data, respectively, and Bern phase. Now all the experimental points, below and above the ρ peak, are inside the allowed range, except the BABAR data below the ρ peak, which are clearly above it. We note that the disagreement is more pronounced for the KLOE 2011 input, and slightly less stringent for the KLOE 2013 input.

V. DISCUSSION AND CONCLUSIONS
In the present paper we have investigated the consistency of the recent high-statistics experimental data on the pion vector form factor with analyticity and unitarity. The aim was to provide additional constraints on the dispersive evaluation of the HVP contribution to the SM value of the muon g − 2. The work was motivated also by the tensions which exist between the most precise data, BABAR and KLOE, which have not been resolved so far.
We have used a model-independent formalism, which uses Fermi-Watson theorem (2) in the elastic region and a very conservative integral condition (3) on the modulus in the inelastic region. From these boundary conditions, rigorous constraints on the values of the form factor at points inside the analyticity domain can be derived, expressed in the particular case of three points as the positivity condition (10) of the determinant (11) and of its minors.
In our analysis, we have exploited this constraint by using as input the experimental data (19) on the form factor on the spacelike axis and the measurements of the modulus in the region (0.65 − 0.71) GeV by the SND, CMD2, BABAR, KLOE 11, KLOE 13 and BESIII experiments. The choice of this region was justified by the fact that data from various experiments are here in better mutual agreement than at higher energies. It may be noted however that this choice is an educated guess and cannot necessarily be considered very rigorous. In particular, some tensions between BABAR and KLOE data exist even in this region.
Using a specific input, we have derived upper and lower bounds on the modulus in the higher-energy region (0.72 − 0.9) GeV, which includes the ρ peak. The bounds derived from analyticity have been converted into central values and standard deviations at 68.2% CL by using numerical simulations on pseudodata samples, and have been presented as a band of "allowed values". The aim was to compare the bands obtained by analytic extrapolation with the experimental data available in the same region, where the tensions between BABAR and KLOE are larger.
We emphasize that in the present formalism the analytic extrapolation is performed without a specific parametrization of the form factor, like the Gounaris-Sakurai formula or the Omnès-like parametrizations used in the literature. Moreover, the results are optimal for a definite input and are rigorously independent on the unknown phase of the form factor above the inelastic threshold. For completeness, we gave in the Appendix A formal argument and a detailed proof of this property, presented here for the first time.
The results presented in Sec. IV indicate some inconsistencies between the data and the extrapolated band. The tension between BABAR and KLOE is manifest and exhibits new facets. As shown in Fig. 4, using as input only BABAR data in the region (0.65 − 0.71) GeV, the KLOE data below the ρ peak are definitely lower than the extrapolated band. More impressively, above the ρ peak all the data, including BABAR, lie below the allowed band. If, on the other hand, only KLOE data in the region (0.65 − 0.71) GeV are used as input, Figs. 6 and 7 show that all the data above 0.72 GeV are consistent with the allowed band, except for BABAR data below the ρ peak, which lie definitely above the band.
A somewhat surprising feature is that all the data above the ρ peak appear to be situated at the lower edge of the extrapolated band. This feature is more dramatic when only BABAR data in the region (0.65 − 0.71) GeV are used as input, as seen in Fig. 4, but is obtained also with input from all experiments, or from all experiments except BABAR, as seen in Figs. 3 and 5, respectively.
It is difficult at present to find a simple explanation for this result. Thus, it is unlikely that the resolution of the puzzle lies in insufficient knowledge of the phase. An important ingredient in our approach is the phase shift δ 1 1 , but this function is calculated with precision in ππ dispersion theory. Moreover, we used two phases, calculated independently in [29] and [30], and the results are consistent, with a somewhat wider bands obtained with Bern phase, which has larger uncertainties near 0.9 GeV.
A global correction factor, with large values precisely in the region of interest, might offer a possible explanation. We recall in this context the so-called ρ − γ mixing factor, calculated in [38], whose effect is to slightly push upward the modulus of the form factor extracted from e + e − data above the ρ peak. However, for the moment we refrain from further speculating on possible explanations of the fact that the data from the region (0.65 − 0.71) GeV seem to require higher values than the experimental measurements of the modulus above the ρ peak, especially above 0.8 GeV.
Instead, we briefly comment on the implications of this result on the HVP contribution to the muon g−2. As it is known, the two-pion LO contribution to a µ , which does not contain the vacuum polarization effects but includes one-photon final-state radiation (FSR), is expressed as an integral over the modulus squared on the pion form factor (for the explicit formulae see for instance Sec. II of [23]).
Using the central value and the error of the extrapolated band shown in Fig. 3, we obtained for the contribution from the region [0.8, 0.9] GeV the value a µ [0.8, 0.9] = (70 ± 4) × 10 −10 . On the other hand, the values obtained for the same quantity from the direct integration of the data are (in units of 10 −10 ) 67.5(4)(6) [5] and 66.6(3) [6], while the fit of the form factor performed in [18] gives 66.6(4). The value obtained by us from analytic extrapolation is much less precise than these estimates, as could be anticipated actually from the figures, and is consistent with them within the uncertainties. However, the fact that the central value is higher by about 3.8 units of 10 −10 is a intriguing and deserves further investigations.
To summarize, in this work we have carried out the logical extension of our previous work which combines the method of unitarity bounds with detailed Monte Carlo simulations to yield determinations of the radius and the values of the form factor below 0.63 GeV, to the region above 0.72 GeV in the present work (as the timelike input data are in the region 0.65−0.71 GeV) which is essentially around the ρ. Since the bounds are now weaker we do not have a high precision determination but rather a detailed consistency test of the compatibility of the data in this region with unitarity and analyticity constraints. Our work reveals some puzzles as discussed in detail in the foregoing. It is a detailed merger of theoretical methods and experimental data in one of the few systems in the low-energy strong interaction sector where high quality data is available in several kinematic regimes.