Pion electromagnetic form factor at high precision with implications to $a_\mu^{\pi\pi}$ and the onset of perturbative QCD

We extend recently developed methods used for determining the electromagnetic charge radius and $a_\mu^{\pi\pi}$ to obtain a determination of the electromagnetic form factor of the pion, $F_\pi^V(t)$, in several significant kinematical regions, using a parametrization-free formalism based on analyticity and unitarity, and with the inclusion of precise inputs from both timelike and spacelike regions. On the unitarity cut, below the first inelastic threshold we use the precisely known phase of the form factor, known from $\pi\pi$ elastic scattering via the Fermi-Watson theorem, and above the inelastic threshold a conservative integral condition on the modulus. We also use as input the experimental values of the modulus at several energies in the elastic region, where the data from $e^+e^-\to \pi^+\pi^-$ and $\tau$ hadronic decays are mutually consistent, as well as the most recent measurements at spacelike momenta. The experimental uncertainties are implemented by Monte Carlo simulations. At spacelike values $Q^2=-t>0$ near the origin, our predictions are consistent and significantly more precise than the recent QCD lattice calculations. The determinations at larger $Q^2$ confirm the late onset of perturbative QCD for exclusive quantities. From the predictions of $|F_\pi^V(t)|^2$ on the timelike axis below 0.63 GeV, we obtain the hadronic vacuum polarization (HPV) contribution to the muon anomaly, $a_\mu^{\pi\pi}|_{\leq 0.63\gev} = (132.97\pm 0.70)\times 10^{-10}$, using input from both $e^+e^-$ annihilation and $\tau$ decay, and $a_\mu^{\pi\pi}|_{\leq 0.63 \gev} = (132.91\pm 0.76)\times 10^{-10}$ using only $e^+e^-$ input. Our determinations can be readily extended to obtain such contributions in any interval of interest lying between $2 m_\pi$ and 0.63 GeV.


I. INTRODUCTION
The electromagnetic form factor of the pion, F V π (t), defined by the matrix element π + (p )|J elm µ |π + (p) = (p + p ) µ F V π (t), where t = q 2 and q = p − p , is a fundamental observable of the strong interactions and a sensitive probe of the composite nature of the pion. An expansion near the origin to linear order in t, F V π (t) = 1 + r 2 π t/6 exhibits the electromagnetic charge radius of the pion, which has recently been determined at high precision in Ref. [1] by a formalism based on analyticity and unitarity with phenomenological input. The result for the electromagnetic charge radius reads r π = r 2 π = (0.657 ± 0.003) fm, which reduced the error by a half from previous determinations. The work was achieved by adapting the techniques introduced in Ref. [2], where the two-pion contribution a ππ µ to the anomalous magnetic moment of the muon was determined in a region where experimental data have significant lack of agreement. In this work, we adapt the methods introduced in these studies to the determination of the form factor itself in several kinematic regions of interest. In contrast to the prior investigations, where a single number was determined in each of them, in the present work we obtain the values of the form factor at a large number of points.
We recall that there is a large amount of information, both theoretical and experimental, on the pion vector form factor, making it one of the most investigated quantity in hadron physics. The form factor determination at high precision is of utmost importance to several observables including the low-energy dipion contribution to the muon g − 2, and poses a significant challenge to experiment as well as to theory. Theoretical studies are based at low energies on nonperturbative approaches and effective theories of the type first formulated by Weinberg [3], and at large energies on perturbative QCD. In the framework of chiral perturbation theory (ChPT), the effective realization of QCD at low energies first formulated at one loop order with two [4] and three light quark flavours [5,6], the pion vector form factor has been calculated up to two loops [7][8][9][10][11]. Lattice gauge theory has recently become another useful nonperturbative tool for the calculation of the form factor at low energies [12,13].
Due to the extensive experimental and theoretical information, the pion vector form factor is, compared with other hadronic quantities, a well-known function. However, the precision does not reach the same level for all timelike and spacelike momenta. A better precision is required on the spacelike axis, for checking the consistency with experimental data and for testing the calculations provided by lattice QCD at low momenta and perturbative QCD at larger momenta. On the timelike axis, at low energies the phase of the form factor is well known, being equal by Fermi-Watson theorem [52,53] to the ππ scattering P -wave phase shift, which has been calculated with high precision using ChPT and Roy equations [54][55][56]. However, the modulus is poorly known, due to the difficulties of the experimental measurements in this region: only two experiments, BABAR [38] and KLOE [40][41][42] reported data at low energies, and unfortunately they are not consistent with each other.
This situation drastically affects the calculation of the hadronic vacuum polarization (HVP) contribution to the muon anomaly, a µ = (g − 2) µ /2, a quantity which plays an important role for testing the standard model and finding possible signals of new physics. The great interest in the muon anomaly is motivated by the present discrepancy of about 3 to 4σ between theory and experiment. New generation measurements of muon g − 2 planned at Fermilab 1 [57] and JPARC [58] are expected to produce results with experimental errors at the level of 16×10 −11 , a factor of 4 smaller compared to the Brookhaven measurement [59]. This requires a precision at the same level also for the theoretical result: see for instance Ref. [58] for an updated review, Refs. [60,61] for most recent phenomenological determinations, and Ref. [62] for a recent lattice calculation.
Dispersion theory, which exploits analyticity and unitarity, is a powerful tool for performing the analytic continuation of the form factor to energies where it is not precisely known. The pion vector form factor is an analytic function in the complex t plane cut along the real axis for t ≥ t + , where t + = 4m 2 π is the first unitarity threshold. Moreover, it is normalized as F V π (0) = 1, and satisfies the Schwarz reflection property F V π (t * ) = (F V π (t)) * . It turns out that the standard dispersion relation, based on the Cauchy integral, is not suitable for F V π (t), since it 1 The E989 experiment at Fermilab has started its pilot runs and is gathering data at an accelerated pace.
requires the knowledge of its imaginary part on the unitarity cut, which is not available in a straightforward way.
On the other hand, as mentioned above, in the limit of exact isospin symmetry, the Fermi-Watson theorem [52,53] states that below the first inelastic threshold, the phase of F V π (t) is equal to the P -wave phase shift of ππ elastic scattering, which is better known. Many dispersion analyses of the pion vector form factor have been based on the so-called Omnès representation, which amounts to reconstruct an analytic function from its phase on the cut. However, this approach involves some assumptions on the phase above the inelastic threshold, where it is not known, and on the positions of the possible zeros in the complex plane. A related approach uses specific parametrizations which implement the analytic properties of the form factor. Recent analyses based on this approach are [63,64].
In the present paper, we use a method based on analyticity and unitarity for calculating the form factor in kinematical regions where it is not precisely known, using the more precise input available in other energy regions. We implement the phase of the form factor along a part of the unitarity cut, where it is well known, and information on the modulus on the remaining part of the cut. Thus, our method is neither a standard dispersive representation, nor a specific parametrization for the analytic extrapolation in the complex momentum plane. The advantage is that we can implement only known input, avoiding to a large extent model-dependent assumptions about the behavior of the form factor in regions where it is less known. The price to be paid was the fact that we do not obtain definite values for the extrapolated quantity, but only optimal allowed ranges for it, in terms of the phenomenological input. This shortcoming has been overcome now as described below.
This method, proposed in [65] and presented in detail in the review [66], has been applied already in several papers [67][68][69][70], where optimal bounds on the pion vector form factor in various energy regions have been derived. An important improvement has been achieved by implementing the statistical distribution of the experimental input by Monte Carlo simulations, which converted the analytic bounds into allowed intervals with definite confidence levels. This elaborate formalism was applied in Refs. [2] and [1] for the calculation with a remarkable accuracy of the low-energy HVP contribution to muon g − 2 and the pion charge radius, respectively. In the present paper, we now complete the task of determining the form factor itself to equally remarkable accuracy.
The outline of this paper is as follows: in Sec. II we review the conditions used as input and formulate the objective of the paper as an extremal problem on a class of analytic functions. In Sec. III we give a detailed description of the information used as input, and in Sec. IV we describe the calculation of the bounds and the Monte Carlo simulation used for implementing the uncertainties of the input data. In this section we also explain the prescription of combining the predictions from different experiments applied in our work. Section V contains our results and Sec. VI a summary and our conclusions. In the Appendix, we present the solution of the functional extremal problem formulated in Sec. II, which is the mathematical basis of our approach.

II. EXTREMAL PROBLEM
Our aim is to make precision predictions for the pion vector form factor in several regions on both spacelike and timelike axis. In particular, we will be interested in the modulus |F V π (t)| in the low energy region t + ≤ t ≤ (0.63 GeV) 2 where t + = 4m 2 π , which will allow a new determination of the pion-pion contribution to the muon anomaly a µ from this region. We will determine also the form factor F V π (t) in the unphysical region 0 ≤ t ≤ t + and at spacelike values t < 0.
We summarize below the conditions adopted as input. We implement first the normalization imposed by gauge invariance at t = 0, expressed by: An important ingredient is Fermi-Watson theorem [52,53] mentioned above. Since this theorem is valid in the exact isospin limit, we must first remove the main isospinviolating effect in the pion vector form factor, known to arise from ω − ρ interference. We shall follow standard approach [71,72] to do this, by defining a purely I = 1 function F (t) as where F ω (t) includes the I = 0 contribution due to ω. Then Fermi-Watson theorem writes as where δ 1 1 (t) is the phase-shift of the P -wave of ππ elastic scattering and t in is the first inelastic threshold.
Above the inelastic threshold t in , where the phase is not known, we shall use the phenomenological information available on the modulus at intermediate energies, and perturbative QCD at high energies. Since the precision is not enough to impose the condition at each t above t in , we shall adopt a weaker condition, written as where ρ(t) > 0 is a suitable positive-definite weight, for which the integral converges and an accurate evaluation of I from the available information is possible. We shall use, in addition, the experimental value of the form factor at one spacelike energy: and the modulus at one energy in the elastic region of the timelike axis, where it is known with precision from experiment: The aim of our work can be expressed as the following functional extremal problem: using as input the conditions (3)-(8), derive optimal upper and lower bounds on |F V π (t)| on the unitarity cut below 0.63 GeV, and on F V π (t) on the real axis for t < t + .
The solution of the extremal problem and the algorithm for obtaining the bounds are presented for completeness in the Appendix. It will be applied in Sec. IV for making precise predictions on the form factor in the regions of interest. In Sec. III we shall describe the phenomenological information used as input.

III. INPUT IN THE EXTREMAL PROBLEM
For the function F ω (t), which accounts for the isospin violation due to ω resonance, we shall use the parametrization 2 proposed in [71,72]: with = 1.9 × 10 −3 . This function is normalized as F ω (0) = 1 and, due to the small value Γ ω = 8.49 MeV [73], is highly peaked around √ t = m ω = 782.65 MeV. In our treatment, we first converted the experimental values of F V π (t) used as input in Eqs. (6), (7) and (8) to the isospin-conserving function F (t) defined in (4), solved the extremal problem for this function and finally reinserted the factor F ω (t) in the results. Actually, since we do not include the resonance region in our study, the corrections due to F ω (t) are very small in all the kinematical regions considered, and are practically negligible for t ≤ 0.
The first significant inelastic threshold t in for the pion form factor is due to the opening of the ωπ channel, i.e., √ t in = m ω + m π = 0.917 GeV. Below this threshold, we use in (5) the phase shift δ 1 1 (t) obtained from dispersion relations and Roy equations applied to ππ scattering in Refs. [54,55] and [56], which we denote as Bern and Madrid phase, respectively. Actually, in the calculation of the Bern phase, for the P -wave phase shift some input from previous data on the form factor was used at the matching point 0.8 GeV, which may raise doubts of a circular calculation (this problem was discussed recently also in [64]). However, we note that the Bern value at 0.8 GeV is practically identical to what has been called "constrained" fit to data (CFD) solution of the Madrid phase [56], which we adopt, and which is independent of form factor data. Actually, the error attached to this input to Bern phase is larger (more than double) than the uncertainty attached to the CFD solution, which reduces the possible bias. Moreover, as we shall explain later, in our determination we take the simple average of the results obtained with the two phase-shifts, which reduces further the potential bias produced by this input and practically avoids the danger of circularity.
We have calculated the integral (6) using the BABAR data [38] 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 1/t decreasing modulus at higher energies, as predicted by perturbative QCD [14,15,20,21]. This choice is expected to overestimate the true value of the integral (see Refs. [67,68,70] for a detailed discussion), which has the effect of leading to weaker bounds due to a monotonicity property discussed in the Appendix. As concerns the weight ρ(t), several choices have been investigated in [70], leading to stable results in most of the investigated regions. In the present work, we have adopted the weight ρ(t) = 1/t, for which the contribution of the range above 3 GeV to the integral (6) is only of 1%. The value of I obtained with this weight is [70] where the uncertainty is due to the BABAR experimental errors. In the calculations we have used as input for I the central value quoted in Eq. (10) increased by the error, which leads to the most conservative bounds due to the monotonicity property mentioned above. On the spacelike axis at moderate and large Q 2 the form factor is extracted indirectly, from experimental measurements of the pion electro-production from a nucleon target, where a virtual photon couples to a pion in the cloud surrounding the nucleon. As a consequence, there are uncertainties associated with the off-shellness of the struck pion and the consequent extrapolation to the physical pion mass pole, which leads to uncertainties in the extraction of the form factor. The errors appear to be under control in the most recent determinations of F π Collaboration at JLab [31,32], as shown in the subsequent analysis [33]. Therefore, as spacelike input (7) we have used the values [31,32] We mention that we do not use as input the data on the spacelike axis near the origin, obtained from eπ scattering by NA7 Collaboration [22]. We shall however compare our predictions for this region with the NA7 data and with the lattice calculations [12].
A major role in increasing the strength of the bounds is played by condition (8). We shall take 0.65 GeV ≤ √ t t ≤ 0.71 GeV, since in this region the modulus measured by various experiments exhibits smaller variations  than in other energy regions and a higher degree of mutual consistency. Moreover, this region is close to the region of interest and therefore has a stronger effect on improving the bounds than the input from higher energies. The e + e − data are taken below 0.705 GeV and the τ -decay data below 0.710 GeV, with the exception of one datum from CLEO that corresponds to an energy of 0.712 GeV. Since this last datum is at an energy that is only marginally higher than the upper limit of the aforementioned energy range, it is included in the analysis.
The numbers of experimental points from various experiments, used as input in our analysis, are summarized in Table I. We emphasize that in this region the e + e −annihilation and τ -decay experiments are fully consistent, so it is reasonable to use all the experiments on an equal footing.
The extraction of the values of timelike modulus |F V π (t)| from the cross-section of the process e + e − → π + π − and the spectral function measured in τ -decay experiments requires the application of several corrections, described in detail in Appendix B of [2]. In particular, for the e + e − experiments the isospin correction due to ω has been applied as discussed above, and the vacuum polarization has been removed from the data.

IV. DETERMINATION OF THE FORM FACTOR AND ITS UNCERTAINTY
Using the algorithm presented in the Appendix, we obtain an allowed range for the value of F V π (t) (or |F V π (t)|) at an arbitrary point t < t in for every set of specific values of the input quantities. However, with the exception of the exact condition F V π (0) = 1, the input quantities are known only with some uncertainties. One of the key aspects of our calculation is the proper statistical treatment of the errors. This is achieved by randomly sampling each of the input quantities with specific distributions: the phase of F V π (t), which is the result of a theoretical calculation, is assumed to be uniformly distributed, while for the spacelike and the timelike data, which are known from experimental measurements, we adopt Gaussian distribution with the measured central value as mean and the quoted error (the biggest error for spacelike data where the errors are asymmetric) as standard deviation. For each point from the input statistical sample, if the input values are compatible, we calculate from Eq. (A16) upper and lower bounds on F V π (t) (or |F V π (t)|). Since all the values between the extreme points are equally allowed, we uniformly generate values of F V π (t) (or |F V π (t)|) in between the bounds. For convenience, the minimal separation between the generated points was set at 10 −3 and for allowed intervals smaller than this limit no intermediate points were created. In this way, for each input from one spacelike t s < 0 and one timelike t t in the region (0.65 − 0.71) GeV, we obtain a large sample of values of F V π (t) (or |F V π (t)|). The results proved to be stable against the variation of the size of the random sample and the minimal separation mentioned above.
In Fig. IV, we present for illustration the distributions of the output values of the form factor at several points of interest (two spacelike points in the upper panel, and two timelike points, one below and the other above the unitarity threshold, in the lower panel). The histograms have been obtained using as input the Bern phase, the value at the spacelike point t s = −1.6 GeV 2 , and the modulus at one timelike point measured by BABAR [38]. Similar results have been obtained using as input the Madrid phase and other experimental data. One can see that the distributions are very close to a Gaussian and allow the extraction of the mean value and the standard deviation (defined as the 68.3% confidence limit (CL) intervals) for the values of interest F V π (t) or its modulus. The next step is to take the average of the results obtained with input from various measurements. Since the degrees of correlations between the measurements at different energies are expected to vary from one experiment to another, we perform first the average of the values obtained with input from each experiment. As argued in [74], the most robust average of a set of n measurements a i is the weighted averagē where δa i is the error of a i . For the best estimation of the error in the case of unknown correlations, the prescription proposed in [74] is to define a function χ 2 (f ) in terms of the covariance matrix C(f ) with elements The parameter f denotes the fraction of the maximum possible correlation: for f = 0 the measurements are treated as uncorrelated, for f = 1 as fully (100%) correlated.
If χ 2 (0) < n − 1, the data might indicate the existence of a positive correlation. The prescription proposed in [74] is to increase f until χ 2 (f ) = n−1. With the solution f of this equation, the standard deviation σ(ā) ofā is determined from the variance [74] On the other hand, a value χ 2 (0) > n − 1 is an indication that the individual errors are underestimated. If the ratio χ 2 (0)/(n − 1) is not very far from 1, the procedure suggested in [73,74] is to rescale the variance σ 2 (ā) calculated with (15) by the factor χ 2 (0)/(n − 1). In our work, this procedure was applied first for combining the results obtained with a definite input phase, a specified input (11) from the spacelike region, and different measurements in the timelike region available from each experiment listed in Table I. In most cases, a large degree of error correlation between the results obtained with different timelike energies was found, as indicated by values close to 1 of the parameter f derived from data. Then the results obtained with the two phases, Bern and Madrid, were combined in a conservative way by taking the simple average of the central values and of the uncertainties. The same conservative average was used for combining the results obtained with the two spacelike data (11).
The last step was to combine the individual values obtained with measurements by the different experiments listed in Table I. Again, the error correlation for these values is difficult to assess a priori. Therefore, we have applied the same data-driven procedure described above for finding the correlations. Since, as discussed in [2], the data from e + e − -annihilation and τ -decay experiments are consistent in the region 0.65 − 0.71 GeV, the results from all the 10 experiments in Table I can be combined into a single central value and standard deviation which we quote as the error.

V. RESULTS
We have applied the procedure described above for deriving central values and standard deviations for F V π (t) in three energy regions: small spacelike momenta Q 2 = −t ≤ 0.25 GeV 2 , where measurements are available from NA7 experiment [22], larger spacelike momenta, up to Q 2 ≤ 8.5 GeV 2 , and the unphysical timelike region 0 < t < t + below the unitarity threshold. We have also derived central values and standard deviations for the modulus |F V π (t)| in the region above the unitarity threshold, below √ t = 0.63 GeV, and have used these results for a new determination of the HVP contribution from energies below 0.63 GeV to the muon g − 2. In the following subsections we present the results for each kinematical region, namely small spacelike momenta, large spacelike momenta, unphysical timelike momenta, and timelike momenta on the unitarity cut below 0.63 GeV. The implications of these determinations are also studied in each of these subsections.

A. Small spacelike momenta
At small spacelike momenta squared, Q 2 ≤ 0.25 GeV 2 , the pion form factor has been measured from ep elastic scattering by the NA7 experiment [22], considered for a long time a landmark experiment. Recently, the ETM collaboration [12] reported the most precise lattice calculations of F V π (−Q 2 ) for small Q 2 . The comparison with the lattice results has been actually the main motivation for choosing this kinematical region in our study. It turns out that our predictions for the form factor in this region are very precise: the errors, obtained by the procedure described in the previous section, vary from 0.0005 near the origin to 0.003 at the end of the region.
In Fig. 2, we present the values of the form factor calculated in this work at a number of spacelike points below 0.25 GeV 2 . Also shown are the experimental data from Ref. [22] and the results of the lattice calculation reported in Ref. [12], shown as a band which includes all the uncertainties. One can see that our results are consistent with the lattice values, and are much more precise. It is a challenge for the future lattice calculations to increase the precision to the level reached by the phenomenological determination based on analyticity and unitarity.
It may be noted that our procedure can be extended further as there is no real constraint on the range of values to be probed in this sector, but for practical purposes, our determination has been limited to the same range as in the lattice study and in the NA7 experiment. The predictions for the pion form factor in the spacelike region near the origin derived in this work, compared with the experimental results of the NA7 Collaboration [22] and the lattice calculations of the ETM Collaboration [12].

B. Large spacelike momenta and the onset of perturbative QCD
It has long been known that in the case of the pion form factor the asymptotic regime described by the dominant term (2) of perturbative QCD sets in quite slowly, due to the complexity of soft, nonperturbative processes in QCD in the intermediate Q 2 region. Several nonperturbative approaches have been proposed for the study of the pion form factor, including QCD sum rules [75], quark-hadron local duality [76][77][78][79], extended vector meson dominance [80], light-cone sum rules [81][82][83], sum rules with nonlocal condensates [84][85][86], AdS/QCD models [87,88], k T factorization method [89], dispersion relations treatment [90], covariant spectator theory [91], and Dyson-Schwinger equation framework [92]. In particular, the onset of the asymptotic regime in the presence of Sudakov corrections [93] and large N c Regge approaches [94] is expected to be quite slow. Constructing a fully valid model to describe the form factor at intermediate energies in fundamental QCD still remains a major theoretical challenge.
Measurements of the spacelike form factor for spacelike momenta are reported in Refs. [23][24][25][26][27][28][29][30][31][32], the most precise being the recent results of the JLab collaboration [31,32] quoted in Eq. (11). The lack of precise experimental data in the higher Q 2 region is a major obstacle to confirm or discard the theoretical models available. The calculation presented in this work provides an alternative way for testing the onset of the asymptotic QCD regime and the validity of various theoretical models proposed for intermediate energies.
In the left panel of Fig. 3, the predictions of this work for Q 2 < 4 GeV 2 , represented as a cyan band which includes the full error, are compared with some of the experimental data. We recall that in our calculation the only input from the spacelike axis consists of the points given in Eq. (11), denoted as Horn in Fig. 3. The increased precision of our determinations is due to the timelike information. One can see that, except for a few points, the experimental measurements are in general agreement with our determinations.
At higher spacelike momenta, the precision of our predictions starts to diminish, since the extrapolation is more sensitive to the values of the form factor at intermediate timelike energies, for which no precise information is available. To account for this, we have adopted the conservative, weaker condition (6). Up to Q 2 around 8 GeV 2 , the precision nevertheless is acceptable, allowing us to probe the onset of the asymptotic regime predicted by factorization and perturbative QCD. In the right panel of Fig. 3, we compare our predictions shown in cyan band with perturbative QCD at LO and NLO, and with some theoretical models. The gray band corresponds to the NLO result obtained by varying the renormalization scale in suitable range following [21] At first sight we note that perturbative QCD at LO can not reliably describe the form factor at Q 2 ≤ 7 GeV 2 . Though the description improves at NLO it is still unreliable for Q 2 ≤ 5.5 GeV 2 . We limit ourselves only to these conservative statements, since precisely at the energies where the NLO and our predictions start to become compatible, our procedure meets its natural limitations. This can be seen in the fact that our band hits the x-axis in right panel of Fig. 3, while there are strong arguments (cf. for instance Ref. [71]) that this form factor cannot have zeros on the spacelike axis. Therefore, we refrain from making definite statements for higher Q 2 , in view of the fact that this is the region where our method lacks the precision that it has in the other three regimes considered in this work.
As we discussed above, there are many theoretical models in the literature for addressing the properties of the form factors in this region. For illustration, we have considered the predictions from four of these as typical examples. For instance, the theoretical models proposed in [78,82] appear to be consistent with the phenomenological band, while the predictions of [88,94] appear to be too high.
We note finally that the results derived in this work are consistent with those derived in our previous work [67], being more precise, since we now included information on the modulus of the form factor on the timelike axis and used extensive Monte Carlo simulations for the error analysis.

C. Unphysical timelike region
No experimental information or QCD lattice calculations are available for the pion form factor in the unphysical timelike region between the origin and the unitarity threshold t + . For this region our method allows to make very precise predictions. In Table II, we list the cen- tral values and the errors on F V π (t) at several unphysical timelike points. This region cannot be accessed by experiment, but it can be by the lattice, in principle, so our results can be viewed as a benchmark for lattice predictions in the future.  In this region the predictions of chiral perturbation theory are expected to be most accurate. The precise determinations in table II can be used to determine the curvature c and higher shape parameter d of the Taylor expansion F V π (t) = 1 + r 2 π t/6 + ct 2 + dt 3 + O(t 4 ).
D. Low energies above the unitarity threshold and the contribution to muon g − 2 As mentioned in the Introduction, above the unitarity threshold, where the form factor is a complex function, its modulus is extracted from the cross section of the e + e − → π + π − process, or, using isospin symmetry, from the hadronic decay of the τ lepton. The τ decay has been for a long time the most precise source of information, in spite of the nontrivial corrections that are required to convert the measured spectral functions into genuine values of |F V π (t)|. However, the accuracy of the e + e − experiments improved gradually, the extraction of the modulus being based at present almost exclusively on them.
At low energies, the modulus of the form factor is poorly known, due to the difficulties of the experimental measurements in this region: only two experiments, BABAR [38] and KLOE [40][41][42] reported data at low energies, and unfortunately they are not consistent among them. Our method allows a precise determination of |F V π (t)| at low energies. In Fig. 4 we present our results, together with the experimental values of BABAR [38] and KLOE [41,42]. For convenience, we show the values of the modulus squared, which enter directly into the calculation of the two-pion contribution to the muon magnetic anomaly. One can see that our predictions are much more precise than the available experimental results, especially at energies below 0.5 GeV. For completeness, we list in Table III the central values and the uncertainties of the modulus squared of the form factor at several energies below 0.63 GeV.  We shall use now these results for making a new determination of the low-energy pion-pion contribution to the muon anomaly. The leading order (LO) two-pion contribution to a µ from energies below √ t up , which does not contain the vacuum polarization but includes onephoton final-state radiation (FSR), is expressed in terms of F V π (t) as In this relation, β 3 π (t) = (1 − 4m π /t) 3/2 is the two-pion phase space relevant for e + e − → π + π − annihilation, is the QED kernel function [95], which exhibits a drastic increase at low t, and is the FSR correction, calculated in scalar QED [96,97].
Using the central values of |F V π (t)| 2 given in Table III, the integral (16) gives (132.97 ± 0.07) × 10 −10 , where we quoted an uncertainty due to the method of integration. In order to estimate the statistical error σ aµ of this result, we shall apply the standard error propagation, expressed in our case as where and Cov(t, t ) is the covariance matrix describing the correlation of the errors on |F V π | 2 at two points t and t . For a most conservative estimate, we assume fully correlated errors, which means that we take where σ(t) is the error on |F V π (t)| 2 , determined by the procedure described in Sec. IV. Then the integral (19) gives an error of 0.69 × 10 −10 . Adding to this the integration error quoted above, we finally obtain a ππ µ | ≤0.63 GeV = (132.97 ± 0.70) × 10 −10 .
For further comparison, we quote also the result obtained using the timelike input on the modulus in the range (0.65-0.71) GeV only from the e + e − experiments: The values (22) and (23) are fully consistent with our previous results (133.26 ± 0.72) × 10 −10 and (133.02 ± 0.77) × 10 −10 , respectively, obtained in [2] for the same quantities with a slightly different method. The difference stems from the fact that in Ref. [2] we generated the statistical distribution of the integral (16) directly from Monte Carlo simulations, without deriving the modulus squared of the form factor at each energy below 0.63 GeV.
We quote also the result a ππ µ | ≤0.63 GeV = 132.5(1.1) × 10 −10 of the recent analysis [64], which exploits analyticity and unitarity by using an extended Omnès representation of the form factor in a global fit of the phenomenological data on e + e − → π + π − from energies below 1 GeV and the NA7 experiment. We note also that the direct integration of the interpolation of the e + e − data below 0.63 GeV, proposed in [61], gives 3 a ππ µ | ≤0.63 GeV = (131.12 ± 1.03) × 10 −10 .
It may be noted that the explicit values listed in Table  III for this region allow an evaluation of the dipion contribution to the muon anomaly in any interval of interest.

VI. DISCUSSION AND CONCLUSIONS
In the present work we have obtained high-precision predictions for the pion electromagnetic form factor in several kinematical regions of interest. We have used a method based on analyticity and unitarity, which does not involve standard dispersion relations or specific parametrizations. The input, summarized in Sec. II, consists of the phase of the form factor on a part of the unitarity cut and a conservative integral condition on the modulus squared on the remaining part. The experimental values at some discrete points on the timelike and the spacelike axes are also included.
Using the solution of a functional extremal problem formulated in Sec. II and discussed in the Appendix, we have derived optimal upper and lower bounds on the values of the form factor (or its modulus) in the regions of interest, which are expressed only in terms of the adopted input and involve no model-dependent assumptions. A key element of our method is the determination of the central values and the errors from statistical distributions obtained from a large set of pseudo-data, and the conservative combination of the results from various experiments using data-driven error correlations. We emphasize that, since we do not use a specific parametrization for the form factor and the analytic bounds do not involve model-dependent assumptions, there are no additional systematic errors in our approach. Thus, the Monte Carlo simulations described in Sec. IV give the full errors of our predictions for the form factor.
We mention that the same technique has been applied in [1] for a precise extraction of the pion electromagnetic charge radius, and in [2] for a direct calculation of the two-pion low-energy contribution to muon g − 2.
The high-precision determinations of the form factor (or its modulus) in several significant kinematical regions are presented in Sec. V. In particular, on the spacelike axis at low Q 2 our results are much more precise than the recent lattice calculations [12], and at larger Q 2 we confirm our previous conclusion [67] that the asymptotic regime of perturbative QCD is away from the region Q 2 ≤ 7 GeV 2 .
On the timelike axis, we derived high-precision values of the modulus squared of the form factor on the unitarity cut below 0.63 GeV, shown in Fig. 4 and Table III. Our predictions are much more precise than the experimental values available in this region from BABAR and KLOE experiments, especially below 0.5 GeV. Also, in Sec. V, the determinations we provide in the unphysical timelike region could serve as a benchmark for theoretical probes in this region.
From the precise values given in Table III, we have performed a new determination of the two-pion contribution from low energies to the muon g − 2. Our results for a ππ µ | ≤0.63 GeV are given in Eqs. (22) and (23), where the first uses the input from both e + e − and τ experiments, and the second only from e + e − experiments. These results are consistent with the values derived in our previous work [2], where the technique of rigorous analytic bounds and Monte Carlo simulations has been applied in a slightly different way, by deriving a statistical distribution directly for the quantity a ππ µ | ≤0.63 GeV . As seen from the values quoted at the end of the previous section, our results are consistent with the prediction of the recent analysis [64] based on analyticity and unitarity, while the result obtained from the direct integration of the data [61] is slightly lower. We emphasize that we do not use as input experimental data from energies below 0.63 GeV or from NA7 experiment. Thus, our prediction for a ππ µ | ≤0.63 GeV is to a certain extent complementary to the determination of the analysis performed in [64], which exploits analyticity and unitarity in a different way and uses as input low-energy data.
This work represents the state of the art in an important low-energy sector of the Standard Model, which is going to be tested at the upcoming Fermilab experiment E969. In contrast to our prior publications [1,2], which were focused on the determination of a single number, here we have obtained an extensive tabulation of the values of the electromagnetic form factor in several significant kinematical regions. Using these determinations, the value of the dipion contribution to the muon anomaly remains consistent with the value reported earlier, proving the robustness of the approach.
The present work combines strong theoretical inputs with modern Monte Carlo methods along with high precision experimental information and phase shift information in regions where experiments are in agreement to shed light on regions where either experiments do not have sufficient precision or where there are significant disagreements, or regions which are not directly accessible by experiment. It also offers a test to theoretical predictions based on very different approaches to strong interaction phenomenology.
From the Fermi-Watson theorem (5), it follows that h(t) is real on the real axis below t in , since the phase of F (t) is exactly compensated by the phase of O(t). Taking into account the fact that h(t) satisfies the Schwarz reflection property, this implies that it is holomorphic on the real axis below t in , having a branch cut only for t ≥ t in .
This relation can be written in a canonical form if we perform the conformal transformation, and express the factors multiplying |h(t)| 2 in terms of an outer function, i.e. a function analytic and without zeros in the unit disk |z| < 1. In practice, it is convenient to construct it as a product of two outer functions [65,66]: the first one, denoted as w(z), has the modulus equal to ρ(t) |dt/dz(t)|. For the choice ρ(t) = 1/t, it is given by the simple expression The second outer function, denoted as ω(z), has the modulus equal to |O(t)|, and can be calculated by the integral representation ) .
This condition leads to rigorous correlations among the values of the analytic function g(z) and its derivatives at points inside the holomorphy domain, |z| < 1 (for a proof and earlier references see Ref. [66]) In particular, in our case this amounts to the positivity condition where the real values z n ∈ (−1, 1) are defined as z n =z(t n ), in terms of the two points t 1 = t s and t 2 = t t used as input and the value t 3 where we want to calculate bounds on the form factor, and ξ n = g(z n ) − g(0), 1 ≤ n ≤ 3.
The inequality (A9) defines an allowed domain for the real values g(z n ). For n = 1 and n = 3 with t < t + , we have from Eqs. (A1) and (A7) while for n = 2 and n = 3 with t > t + g(z n ) = w(z n ) ω(z n ) |F (t n )|/|O(t n )|, where the modulus |O(t)| of the Omnès function is obtained from (A2) by the principal value (PV) Cauchy integral . (A15) One can show that for each specific input, the determinant (A10) is a concave quadratic function of the unknown value F (t 3 ) for t 3 < t + , or the modulus |F (t 3 )| for t 3 > t + , so the inequality (A9) can be written as where x = F (t 3 ) or x = |F (t 3 )|. This inequality leads to a definite allowed range for x if B 2 − AC ≥ 0 and has no solution if B 2 − AC < 0. The latter case occurs when the phenomenological input adopted is inconsistent with analyticity. From the inequality (A16), one can obtain upper and lower bounds on F (t 3 ) (or |F (t 3 )|), expressed in terms of the adopted input. Finally, the isospin correction is applied back according to (4), for obtaining the desired bounds on the form factor F V π (t). One can prove [65,66], that the bounds are optimal and their values do not depend on the unknown phase of the form factor above the inelastic threshold t in (the dependence of the Omnès function (A2) on the arbitrary phase δ(t) for t > t in is compensated exactly by the corresponding dependence of the outer function (A6)). Furthermore, for a fixed weight ρ(t) in (6), the bounds become stronger/weaker when the value of the value of I is decreased or increased, respectively. These properties make the formalism model independent and robust against the uncertainties from the high energy region.