Dispersive analysis of the experimental data on the electromagnetic form factor of charged pions at spacelike momenta

The experimental data on the electromagnetic form factor of charged pions available at spacelike momenta are analyzed using the Dispersive Matrix (DM) approach, which describes the momentum dependence of hadronic form factors without introducing any explicit parameterization and includes properly the constraints coming from unitarity and analyticity. The unitary bound is evaluated nonperturbatively making use of the results of lattice QCD simulations of suitable two-point correlation functions contributing to the HVP term of the muon. Thanks to the DM method we determine the pion charge radius from existing spacelike data in a completely model-independent way and consistently with the unitary bound, obtaining $_{DM} = 0.703 \pm 0.027$ fm. This finding differs by $\simeq 1.6$ standard deviations from the latest PDG value $_{PDG} = 0.659 \pm 0.004$ fm, which is dominated by the very precise results of dispersive analyses of timelike data coming from measurements of the cross section of the $e^+ e^- \to \pi^+ \pi^-$ process. We have analyzed the spacelike data using also traditional $z$-expansions, like the Boyd-Grinstein-Lebed (BGL) or Bourrely-Caprini-Lellouch (BCL) fitting functions and adopting a simple procedure that incorporates ab initio the non-perturbative unitary bound in the fitting process. We get $_{BGL} = 0.711 \pm 0.039$ fm and $_{BCL} = 0.709 \pm 0.028$ fm in nice agreement with the DM result. We have addressed also the issue of the onset of perturbative QCD by performing a sensitivity study of the pion form factor at large spacelike momenta, based only on experimental spacelike data and unitarity. Hence, although the leading pQCD behaviour is found to set in only at very large momenta, our DM bands may provide information about the pre-asymptotic effects related to the scale dependence of the pion distribution amplitude.


I. INTRODUCTION
Since the pion is the lightest bound state in QCD, its physical properties carry important information about the way quark and gluon degrees of freedom govern the low-energy dynamics.Therefore, its precise determination represents an important test for our fundamental theory of the strong interactions and requires nonperturbative theoretical approaches, like QCD simulations on the lattice.The electromagnetic (em) form factor of a (charged) pion, F V π (Q 2 ), is defined in pure QCD by the matrix element where q = p − p ′ is the 4-momentum transfer, Q 2 ≡ −q 2 and J em µ is the em current operator, namely with q f being the electric charge of the quark with flavor f in units of the electron charge.
For spacelike values of the squared 4-momentum transfer (Q 2 ≥ 0 or, equivalently, q 2 ≤ 0) the em pion form factor contains information on the distribution of its charged constituents, namely valence and sea light quarks, while for timelike values it has a branch cut starting at the annihilation threshold 4M 2 π .For Q 2 ≤ −4M 2 π (q 2 ≥ 4M 2 π ) it becomes complex and its modulus is a crucial quantity governing the 2π contribution to the hadronic vacuum polarization (HVP) of the muon anomalous magnetic moment (see, e.g., Ref. [3]).As well known, the muon HVP has long played an important role for testing the Standard Model of particle physics.
The experimental information on the em pion form factor is quite rich.At spacelike values of Q 2 the form factor has been determined using electron-pion scattering experiments [4][5][6] and pion production off nucleons [7][8][9][10][11][12][13][14].In the timelike region the modulus of the pion form factor has been extensively measured using the cross section of the process e + e − → π + π − (see Ref. [3] for a recent compilation) as well as data on the hadronic τ decays in the limit of isospin symmetry.Concerning the extraction of the em pion form factor from experimental data and its analysis in terms of dispersion methods a consistent treatment of radiative corrections (due to both vacuum polarization and final-state radiation effects) must be guaranteed, as described in Refs.[15,16].In the case of the spacelike data the radiative corrections considered in the experiments include already the subtraction of vacuum polarization effects [17,18].Thus, the dispersive treatment can be applied to the spacelike data for the em pion form factor without any adjustment.
An important quantity characterizing the em pion form factor is its slope at Q 2 = 0, more precisely the pion charge radius, ⟨r π ⟩, defined as ⟨r π ⟩ ≡ ⟨r 2 π ⟩ with ⟨r 2 π ⟩ ≡ −6 where the rightmost formula can be obtained via dispersion relations.In the latest PDG review [2] the result for ⟨r π ⟩ reads ⟨r π ⟩ P DG = 0.659 ± 0.004 fm , coming from an average of four different results: ⟨r π ⟩ = 0.656 ± 0.005 fm representing a suitable average of the analyses of timelike (e + e − ) and spacelike [5] data made in Refs.[16,19], ⟨r π ⟩ = 0.663 ± 0.023 fm using the spacelike data from the F2 experiment at FNAL [4], ⟨r π ⟩ = 0.663 ± 0.006 fm using the spacelike data from the NA7 experiment at CERN [5] and ⟨r π ⟩ = 0.65 ± 0.08 fm using the spacelike data from the SELEX experiment at FNAL [6].The first result from Refs.[16,19] is based on a dispersive representation of the em pion form factor, which properly satisfies unitarity and analyticity.On the contrary, the other three results, based only on spacelike data, are obtained by fitting the data with a simple monopole Ansatz, which may introduce a disturbing model dependence and may be inconsistent with unitarity (see similar remarks made in Ref. [20], where Padé approximants are employed).
The aim of the present work is to describe the Q 2 -dependence of the experimental data on the em pion form factor at spacelike momenta without introducing any explicit parameterization and fulfilling at the same time the constraints coming from unitarity and analyticity.This will allow us to determine the pion charge radius from existing spacelike data in a completely model-independent way, while fulfilling unitarity.This goal can be achieved by adopting the Dispersive Matrix (DM) approach developed in Ref. [1], and already applied successfully to the description of the hadronic form factors relevant in semileptonic B-meson weak decays in Refs.[21][22][23][24][25].In this work the unitary bound on F V π (Q 2 ) will be imposed using for the first time a nonperturbative determination of the relevant transverse vector susceptibility obtained using the results of lattice QCD (LQCD) simulations of suitable two-point correlation function contributing also to the muon HVP.Our result is ⟨r π ⟩ DM = 0.703 ± 0.027 fm , which differs by ≃ 1.6 standard deviations from the PDG value (4) with an uncertainty much larger (by a factor ≃ 4.5) than the one quoted in the experimental work of Ref. [5].We have analyzed the spacelike data using also traditional z-expansions, like the renowned Boyd-Grinstein-Lebed (BGL) [26] or Bourrely-Caprini-Lellouch (BCL) [27] fitting functions and adopting a simple procedure that easily incorporates ab initio the non-perturbative unitary bound in the fitting process.We get ⟨r π ⟩ BGL = 0.711 ± 0.039 fm and ⟨r π ⟩ BCL = 0.709 ± 0.028 fm in nice agreement with the DM result (5).The analysis of the em pion form factor using the basic features of the DM approach is not new at all.In Refs.[19,[28][29][30][31] the authors adopted a matrix approach similar to the DM one, characterized by the use of one timelike and one spacelike constraint at the same time, and by a subsequent suitable averaging procedure of the results corresponding to different pairs of input data.In this work we introduce a new procedure, the unitary sampling procedure, valid for any number of data points.In this way we demonstrate that the DM approach is an easy and very effective tool for analyzing even large sets of data points (more than 50 in this work) fulfilling exactly the unitarity and analyticity constraints.
The structure of this work is as follows.
In Section II and Appendix A we summarize the main features of the DM approach applied to the em pion form factor.In particular, we elucidate the meaning of the DM unitary filter, which allows to select in a model-independent way only the subset of input data that can be reproduced exactly by a unitary z-expansion.This feature is not guaranteed by approaches based on explicit z-expansions and it becomes more important as the impact of the unitary filter is more severe.In these approaches the attention is focused only on the fitting function and not also on the fitted data (either experimental or theoretical ones).Even if the fitting function is constructed to satisfy unitarity, the fitting procedure is applied to all input data regardless whether the latter ones satisfy unitarity or not (i.e., regardless whether the input data can be exactly reproduced by a unitary z-expansion).We point out that fitting non-unitary input data might introduce distortions in unitary z-expansions (see Section VII for a numerical evidence).Up to our knowledge this potential problem is avoided only in the DM method.
In Section III we discuss the non-perturbative determination of the unitary bound used in this work, namely the transverse vector susceptibility χ T , obtained using the results of lattice QCD (LQCD) simulations of suitable two-point correlation functions contributing to the HVP term of the muon.
In Sections IV and V we apply the DM method to the electroproduction JLAB-π data [13].Since the unitary bound turns out to be extremely selective as the number of data points increases, we develop an efficient procedure to generate a distribution of values for the pion form factor satisfying unitarity, i.e. to get a set of unitary input data, valid for any number of data points.The unitary sampling procedure is described in detail in the case of the electroproduction JLAB-π data and it can be easily generalized to any set of hadronic form factors, which must satisfy unitary bounds.
In Section VI the unitary sampling method is applied to both the CERN [5] and electroproduction JLAB-π [13] data for a total of more than 50 data points.The DM band for the em pion form factor is positively compared with the results obtained in Ref. [16] by means of a unitary analysis of both timelike e + e − and spacelike CERN data.A difference is observed at small values of Q 2 , which translates into the value (5) of the pion charge radius w.r.t. the result ⟨r π ⟩ = 0.655 ± 0.003 fm from Ref. [16].
In Sections VII and VIII we analyze the spacelike data using the BGL [26] and BCL [27] z-expansions, respectively.We adopt a simple procedure that incorporates ab initio the non-perturbative unitary bound, described in Appendices B and C. A detailed comparison among the unitary BGL fitting procedure and the DM method is performed, showing explicitly that distortions are produced at large spacelike values of Q 2 by fitting non-unitary input data.
In Section IX we investigate the role of the auxiliary quantity Q 2 0 , at which the transverse susceptibility χ T (Q 2 0 ) is evaluated, on the unitary DM filter and on the corresponding DM band for the em pion form factor.
In Section X we address the issue of the onset of perturbative QCD (pQCD) at large spacelike values of Q 2 .As well known, for Q 2 → ∞ the leading behaviour of F V π (Q 2 ) predicted by pQCD [32][33][34][35] is given by 8πf 2 π α s (Q 2 )/Q 2 , where f π ≃ 130 MeV is the pion decay constant and α s (Q 2 ) is the running strong coupling.We perform a sensitivity study and present the DM predictions for Q 2 ≳ 5 GeV 2 based only on unitarity and experimental data available at spacelike momenta.Although the leading pQCD behaviour is found to set in only at very large momenta, our DM bands may provide information about the pre-asymptotic effects related to the scale dependence of the pion distribution amplitude.
Our conclusions are summarized in Section XI.We point out that the DM approach is equally well suited to be applied also to available results of LQCD calculations of F V π (Q 2 ) and, more generally, to experimental plus LQCD data on F V π (Q 2 ).In this work, since tensions are present among e + e − experiments (further exacerbated by the recent results from the CMD-3 Collaboration [36]), we are interested in the analysis of the experimental data available at spacelike momenta without any mixing with timelike data, allowing in this way an interesting comparison with the results of Ref. [16], which are based almost totally on timelike data.We leave the DM analysis of LQCD data as well as of timelike plus spacelike data to future separate works.

II. THE DM APPROACH FOR THE EM PION FORM FACTOR
The DM approach is a non-perturbative method for computing hadronic form factors in a model-independent way in their full kinematical range [1,37].
The starting point is a dispersive bound that, for a generic form factor f , can be written as [26,38,39] where ϕ(z) is a kinematical function dependent on the specific spin-parity channel and χ is the so-called susceptibility, related to the derivative of the Fourier transform of a suitable Green function of bilinear quark operators [26].The conformal variable z(t) is defined as where t = q 2 = −Q 2 is the squared 4-momentum transfer and, for the case of interest in this work, t + = 4M 2 π and t 0 = 01 .
In the case of sub-threshold bound-state poles located at t i = M 2 Ri < t + , the requirement of analyticity can be fulfilled by modifying the kinematical function ϕ(z) through the so-called Blaschke factors, namely [37] , where z(t) is the complex conjugate of the conformal variable z(t).In the case of the em pion form factor no sub-threshold pole is present.By introducing the inner product [37,40] where ḡ(z) is the complex conjugate of the function g(z), Eq. ( 6) can be also written as Following Refs.[37,40] we introduce the set of functions so that the use of Cauchy's theorem yields .
The central ingredient of the DM method is the matrix [37,40] where t 1 , . . ., t N are the values of the squared 4-momentum transfer at which the form factor f (z) is known.Note that the DM method can be applied not only to a series of theoretical values f (z(t i )) (with i = 1, 2, ...N ), but also directly to experimental data (as done in this work and in Ref. [41]).
The important feature of the matrix M is that, thanks to the positivity of the inner products, its determinant is positive semidefinite, i.e. det M ≥ 0. This property is not modified when the matrix element ⟨ϕf |ϕf ⟩ is replaced by the upper bound given by the susceptibility χ through Eq. ( 8).Thus, the original matrix (9) can be replaced by where ϕ i f i ≡ ϕ(z i )f (z i ) (with i = 1, 2, ...N ) represent the known values of ϕ(z)f (z) corresponding to the given set of values z i .By imposing the positivity of the determinant of the matrix (10) it is possible to explicitly compute the lower and upper bounds that unitarity imposes on the form factor f (z) for a generic value of z on the real axis, namely [1] where When z → z i one has d(z) ∝ 1/(z − z i ) and, therefore, β(z) → f i and γ(z) → 0. In other words, Eq. ( 11) exactly reproduces the set of input data {f i }.In a frequentist language this corresponds to a vanishing value of the χ 2 -variable.
Unitarity is satisfied only when γ(z) ≥ 0, which implies the condition χ ≥ χ DM .Such a condition depends on the set of input data {f i } and it is independent on any parameterization or fitting Ansatz of the input data.
The meaning of the DM filter χ ≥ χ DM is clearer in terms of explicit z-expansions, like the BGL ones [26].When χ ≥ χ DM , it is guaranteed the existence of (at least) one BGL fit (either truncated or untruncated) that satisfies unitarity and, at the same time, reproduces exactly the input data.On the contrary, when χ < χ DM , a unitary z-expansion passing through the data does not exist, since the input data do not satisfy unitarity.The important feature of the DM approach is that only the unitary input data are eligible for consideration, while those data that do not satisfy the unitary filter χ ≥ χ DM are discarded.
We want to elucidate better the relevance of the above feature of the DM approach.Let us consider a sample of input data corresponding to known values of the form factor at a series of points z i generated according to a given covariance matrix.For each event we can apply the DM filter χ ≥ χ DM and, consequently, we can divide the original sample into two disjoint subsets: the one corresponding to input data satisfying the DM filter and the one made of non-unitary events.In what follows we will refer to the first subset as the unitary input data and to the second one as the non-unitary input data.
The DM approach provides a band of values for the form factor f (z) which are consistent with unitarity without making use of any explicit fitting Ansatz.More precisely, the DM band is given by the convolution of the uniform distribution corresponding to Eqs. ( 11)- (13) with the distribution of the unitary input data {f j }, i.e.only those fulfilling the condition χ ≥ χ DM (see later Section V).In other words, the DM approach automatically provides the envelope of the results of all possible (either truncated or untruncated) z-expansions, which satisfy unitarity and at the same time exactly reproduce the unitary input data.
We stress again that separating the input data in the two disjoint subsets corresponding to either χ ≥ χ DM or χ < χ DM is an important feature of the DM approach, which is not guaranteed by approaches based on explicit z-expansions (including the one of Ref. [42] and also the recent Bayesian approach of Ref. [43]).Indeed, in these approaches the attention is focused only on the fitting function and not also on the fitted data (either experimental or theoretical ones).Even if the fitting function is constructed to satisfy unitarity, the fitting procedure is applied to all the input data regardless whether they satisfy unitarity or not (i.e., regardless whether the input data can be exactly reproduced by a unitary z-expansion).In the case of the unitary subset of input data it is always possible to find a suitable BGL fit, that satisfies unitarity and at the same time exactly reproduces the input data.This corresponds to the possibility to reach a null value of the χ 2 -variable by increasing the order of the truncation of the BGL fit (up to the number of data points).On the contrary, when the input data do not satisfy the unitary filter, it is not possible to find a fitting z-expansion that satisfies unitarity and at the same time exactly reproduces the input data.This corresponds to a non-vanishing value of the χ 2 -variable, which depends on the impact of the non-unitary input data.The above considerations applies equally well also to the case of explicit z-expansions like the BCL ones [27].
It is clear that the application of a fitting function (even if unitary) to a subset of input data that do not satisfy unitarity may lead to a distortion of the fitting results related directly to the impact of the non-unitary effects present in the input data.In particular, such a distortion may be relevant when the fitting function extrapolates the form factor in a kinematical region not covered by the input data.Thus, the application of the DM method is simpler and more general w.r.t.other approaches, like the BGL or BCL z-expansions, particularly when the number of input data increases and the unitary constraint becomes more selective.In these cases a BGL or BCL truncated expansion would require to take into account a large number of fitting parameters with no guarantee of avoiding the non-unitary effects possibly present in the input data.This issue will become evident in Sections VII and VIII, where we apply unitary BGL or BCL approaches to analyze the spacelike data for the pion form factor.In Section VII we address explicitly the issue of the impact of non-unitary input data on a unitary BGL fit.
Let us now consider explicitly the case of the em pion form factor assuming that it is known for a series of (N + 1) values Q 2 i , namely where we have added the value i = 0 to include the absolute normalization condition F V π (Q 2 0 = 0) = 1.The susceptibility relevant for the em pion form factor is the one of the transverse vector channel χ T (see next Section for its explicit definition) and we introduce an auxiliary variable Q 2 0 at which the susceptibility χ T is evaluated.Denoting by z i the value of the conformal z-variable corresponding to Q 2 i , i.e.
π , the constraint due to unitarity and analyticity on the values F i can be written in the form where and the kinematical factor ϕ i (Q ) is explicitly given by [26,44] with For a generic value of z on the real axis, when the unitary filter (18) is satisfied, the pion form factor F V π (z) is limited by the bounds where ϕ(z, Q while χ DM (Q 2 0 ) and z 0 are given by Eqs. ( 19) and (22), respectively.We stress once more that the important feature of the DM method is the possibility to predict the value of the form factor F V π (Q 2 ) at a generic value of Q 2 using only the knowledge of the pion form factor at the series of values Q 2 i without any reference to a specific parameterization, provided the unitary filter ( 18) is fulfilled.We will describe the DM procedure in detail later on in Section V, while in the next Section we address the non-perturbative determination of the transverse vector susceptibility χ T (Q In QCD the transverse vector susceptibility χ T (Q where V 2π (τ ) is the 2π contribution to the Euclidean vector-vector current correlator V (τ ), τ is the Euclidean time distance and j 1 (x) is an ordinary Bessel spherical function.Note that at Q 2 0 = 0 the susceptibility χ T (Q 2 0 = 0) is proportional to the fourth moment of the correlator V 2π (τ ), which contributes to the Hadronic Vacuum Polarization term of the muon (g − 2) (see Ref. [3]).
As well known, the Euclidean correlator V (τ ) can be obtained by taking the Fourier transform of the spatial components of the HVP tensor, which in turn is related via dispersion relations to the (one photon) e + e − annihilation cross section into hadrons, namely (see, e.g., Ref. [45]) where with ω being the center-of-mass energy.
In QCD, neglecting the electron mass, the 2π contribution R 2π (ω) to R had (ω) is given by (see, e.g., Ref. [16]) where F V π (ω) is the em pion form factor in the time-like region q 2 = ω 2 ≥ 4M 2 π .The 2π contribution to the Euclidean correlator reads as Since the transverse vector susceptibility χ T (Q Note that: • the integrand in the r.h.s. of Eq. ( 34) is positive definite at all energies; • for large values of Q 2 0 the transverse susceptibility drops down as fast as 1/Q • Eq. ( 34) can be cast in the form of the unitary bound (6) as a strict equality.Indeed, the relation between the conformal variable z on the unit circle |z| = 1 and the center-of-mass energy ω is so that one has ω = 2M π 2/(1 − cosα).It follows that where ϕ(z, Q 2 0 ) is the kinematical function given in Eq. ( 27).The Euclidean correlator of two em currents has been evaluated on the lattice by several collaborations (see, e.g., Ref. [3]).In particular, the 2π contribution V 2π (τ ) has been estimated at the physical point and in the continuum and infinite volume limits in Ref. [46] 2 .Thus, using the lattice-based correlator V 2π (τ ) we have evaluated the quantity 4M 2 π χ T (Q 2 0 ) using Eq. ( 28).The Q 2 0 dependence obtained in this way is shown in Fig. 1 by the blue dots.Alternatively, we have calculated Eq. ( 34) adopting for |F V π (ω)| the results of the dispersive analysis of the e + e − data available from Ref. [16] up to ω = 1 GeV and putting |F V π (ω)| = 0 for ω > 1 GeV.We will refer to the results obtained in this way for the susceptibility 4M 2 π χ T (Q 2 0 ) as the data-driven ones, represented in Fig. 1 by the black line.Since the pion form factor of Ref. [16] is provided up to ω = 1 GeV and the integrand in the r.h.s. of Eq. ( 34) is positive definite, the data-driven results for 4M 2 π χ T (Q 2 0 ) represent a lower bound for the transverse susceptibility.Reassuringly, the lattice-based results turn out to be slightly higher than the data-driven ones.The agreement is remarkable with differences not exceeding ∼ 10% up to The agreement shown in Fig. 1 can be improved by adding to the data-driven pion form factor a simple power-law tail for ω > 1 GeV of the form 4 , which reproduces within the errors the results of Ref. [16] in the ω 2 -range from ≃ 0.9 GeV 2 up to 1 GeV 2 .In this case we have found that the differences w.r.t. the lattice-based results for the susceptibility 4M 2 π χ T (Q 2 0 ) do not exceed ≃ 4% (i.e., less than ≃ 2 standard deviations).It should be kept in mind, however, that the parameterization of the pion form factor F V π (ω) adopted in Ref. [16], while fulfilling the requirements of unitarity and analyticity, includes the contributions from 2π, 3π and inelastic ) versus Q 2 0 evaluated using Eq. ( 28) and adopting the 2π correlator V2π(τ ) obtained in Ref. [46] using LQCD simulations (see text).Black line: the susceptibility 4M 2 π χT (Q 2 0 ) evaluated using Eq. ( 34) adopting for |F V π (ω)| the results of the dispersive analysis of the e + e − data available from Ref. [16] up to ω = 1 GeV and putting |F V π (ω)| = 0 for ω > 1 GeV.The uncertainties of the data-driven results are hardly visible on the vertical scale being of the order of ≃ 0.6%. .channels, estimated conservatively up to ω = 1 GeV.Any extension to the ω-region above 1 GeV requires at least the inclusion of the contributions arising from higher resonances like ρ(1450) and ρ(1700) (see the corresponding note in the PDG review [2]), which is still to be settled [16].At the same time, also the application of the approach of Ref. [46] for estimating the 2π correlator V 2π (τ ) from lattice isoQCD simulations is limited to Euclidean time distances τ ≳ 1 fm, which qualitatively corresponds to energies below ≈ 1 GeV (see Eq. ( 32)).According to Eq. ( 34) the impact of the high-energy tail of F V π (ω) increases as Q 2 0 increases and, therefore, we consider trustable the results obtained for both the lattice-based and the data-driven transverse susceptibility only up to Q 2 0 ≈ 1 GeV 2 (as adopted in Fig. 1).In what follows, we will make use of the lattice-based results for the transverse vector susceptibility 4M 2 π χ T (Q The main reason is that we want to analyze the spacelike data within the DM method without using information coming from e + e − data.In this way we will compare our results coming only from the spacelike sector with the corresponding ones obtained in Ref. [16] from timelike data (see later Section VI). At This value can be compared with the upper bound provided by the quantity 4M , where Π is the isovector contribution to the slope of the HVP polarization function evaluated at vanishing four-momentum transfer.The isovector HVP slope, which contains contributions also from intermediate states other than the 2π states, has been calculated (in isoQCD) by several lattice collaborations, namely BMW [55], RBC [56] and FHM [57], obtaining, respectively, 4M 2 π Π (I=1) 1 = 0.00607 (19), 0.00624 (17), 0.00611 (9).

IV. ANALYSIS OF THE ELECTROPRODUCTION DATA
Presently the experimental data on the em pion form factor at spacelike momenta can be divided into two groups.For values of Q 2 ≲ 0.25 GeV 2 the pion form factor F V π (Q 2 ) has been determined by measuring the scattering of high-energy (on-shell) pions off atomic electrons at FNAL F2 [4], CERN SPS [5] and FNAL SELEX [6].The data are shown in the upper panel of Fig. 2.
At higher values of Q 2 the pion form factor is extracted from cross section measurements of the reaction 1 H(e, e ′ π + )n, i.e. from pion electroproduction off the proton, which implies initial off-shell pions.In such a case Experimental data on the em pion form factor F V π (Q 2 ) obtained at FNAL F2 [4] (red upper triangles), CERN SPS [5] (black circles), FNAL SELEX [6] (green lower triangles) and from the reanalysis of the electroproduction data performed by the JLAB-π Collaboration in Ref. [13] (blue squares).
the separation of the longitudinal and transverse response functions as well as the extrapolation of the observed scattering from virtual pions to the one corresponding to on-shell pions have to be carefully considered for estimating the systematic uncertainties.Using the electroproduction technique the pion form factor F π V (Q 2 ) has been determined for Q 2 ≳ 0.35 GeV 2 in various experiments at CEA/Cornell [7], DESY [8,9] and JLAB [10][11][12]14].
A careful analysis of the systematic uncertainties of all the electroproduction data was carried out by the JLAB-π Collaboration in Ref. [13].Particular attention was paid to estimating the uncertainty due to the extrapolation of the electroproduction data to the pion pole.Within the errors no inconsistency is visible with the pion form factor obtained by the dispersive analysis of the e + e − data of Ref. [16] (see later on Fig. 10) as well as with available results from lattice (iso)QCD simulations (see the review in Ref. [58]).
The results of the JLAB-π collaboration are collected in Table I together with the absolute normalization value F 0 = 1 at Q 2 0 = 0 (i.e.z 0 = 0) corresponding to the charge conservation.The electroproduction data are shown in the lower panel of Fig. 2 together with the CERN ones.Note that: i) few experimental data are plagued by large and asymmetric errors, which come from the systematic uncertainty due to a long extrapolation to the pion pole [13]; ii) at Q 2 = 1.6 GeV 2 Ref. [13] quotes two results: 0.233 +19 −17 and 0.243 +23 −14 .Both are shown in the lower panel of Fig. 2, while in Table I only their average 0.238 +21 −17 is considered.We now apply the DM approach to the set of N + 1 = 8 data collected in Table I assigning a very small, but  I. Experimental data on the em pion form factor F V π (Q 2 ) obtained using the electroproduction technique analysed by the JLAB-π Collaboration in Ref. [13] together with the absolute normalization value F0 = 1 at Q 2 0 = 0 (i.e.z0 = 0) due to charge conservation.The fifth column represents the data after symmetrization of the errors.The sixth column contains the values of the kinematical coefficients di (see Eq. ( 20)).The last column shows the values of the kinematical factors ϕi evaluated at Q 2 0 = 0 (see Eq. ( 21)).
non-vanishing error to the data point F 0 = 1 at Q 2 = 0, namely σ 0 = 10 −16 .We also symmetrize the errors obtaining the set of data points shown in the fifth column of Table I.Since no information is available on the covariance matrix of the electroproduction data, the form factor values are considered to be uncorrelated, namely the covariance matrix C is given by where σ 2 i is the variance of F i with i = 0, 1, ...N .
We start by choosing Q 2 0 = 0 and postponing to Section IX the discussion about the impact of a generic choice Q 2 0 > 0. We assume a gaussian distribution for the non-perturbative transverse susceptibility 4M 2 π χ T (Q 2 0 = 0) = 0.00574 (10) (see Eq. ( 38)).This distribution is taken to be uncorrelated with those of the form factor points collected in Table I.
A sample of 10 5 uncorrelated events normally distributed is generated using as input the mean values F ≡ {F i } and uncertainties σ ≡ {σ i } with i = 0, 1, ...N .For each event we calculate the susceptibility χ DM (Q 2 0 = 0) given by Eq. (19).It turns out that the calculated values of χ DM (Q 2 0 = 0) range from a minimum equal to ∼ 2.6 up to a maximum given by ∼ 1.4 • 10 9 and, therefore, none of the 10 5 generated events satisfies the unitary filter (18).The same happens also when we increase the size of the sample up to 10 6 .
The main reason for the above finding can be traced back to the values of the kinematical coefficients d i , given by Eq. (20).These coefficients depend only on the series of values z i and their numerical values are shown in the sixth column of Table I.They turn out to be quite large in absolute value and to have alternating signs.It is therefore very unlikely to generate an event with uncorrelated values of the form factor points leading to a value of χ DM (Q A very delicate compensation among the contributions of the various data points to Eq. ( 19) is required and this naturally implies specific correlations among the form factor points.In principle, one may increase the size of the sample until some of the events satisfy the unitary filter, but a brute-force increase of the size of the sample may become impracticable for large values of the number of data points N (see also later on Section VI).
As already pointed out in Section II, the unitary filter (18) acts as a constraint and it allows to select a subset of the initial events F and C ≡ {C ij }, made only by unitary events.Such a subset corresponds to new values F for the form factor points and to a new covariance matrix C representing the unitary form factor points (and correlations) on which any further analysis fulfilling unitarity must be based.Thus, we have to find an efficient way to determine the unitary values F and C. In the next Section we illustrate a simple procedure able to achieve this goal and applicable for any value of N .

V. UNITARY SAMPLING PROCEDURE
The gaussian multivariate distribution used in the previous Section is based on the probability density function (PDF) given by where {F i , } and {C ij } are respectively the mean values and the covariance matrix used as inputs.As well known, the PDF (40) favors the relative likelihood of small values of the quadratic form , which however may correspond to large values of the susceptibility (19), as shown in the previous Section.
We now modify the above PDF in order to allow the susceptibility (19) to be small enough to fulfill the unitary constraint (18).We consider the following new PDF: where s is a parameter (expected to determine the number of events satisfying the unitary filter ( 18)) and the matrix ) is defined as The use of Eq. ( 41) as a PDF allows to increase the relative likelihood of small values of the susceptibility χ DM (Q 2 0 ) at the expense of decreasing the PDF (40).Introducing the matrix C defined in compact notation as Eq. ( 41) can be easily rewritten in the form where the new vector of mean values F is related to the starting one F by Note that the second exponential in the r.h.s. of Eq. ( 44) does not depend on {f i } and therefore it is irrelevant for the relative likelihood of the events, so that the new PDF is simply given by which represents a multivariate gaussian distribution characterized by the new set of input values { F i , C ij } given by Eqs. ( 45) and (43), respectively.In the case of the electroproduction data of Table I we generate samples with 10 5 events according to the new PDF (46) for various values of the parameter s (assuming again Q 2 0 = 0).Then, we calculate the corresponding susceptibility χ DM (Q 2 0 = 0) given by Eq. (19).The results are shown in Fig. 3, where the (gaussian) distribution corresponding to the non-perturbative transverse result 4M  19)) obtained using samples of 10 5 events generated according to the multivariate distribution (46) for various values of the parameter s.The quantity 0.00574 represents the central value of the non-perturbative transverse susceptibility (38), obtained in Section III at Q 2 0 = 0.The percentage p of events passing the unitary filter ( 18) is given in the inset for each value of s.The grey histogram represents the gaussian distribution of the non-perturbative transverse susceptibility ratio 4M 2 π χT (0)/0.00574.
be clearly seen that, as the parameter s increases, the PDF ( 46) can be very effective in generating events with χ DM (Q satisfying the unitary filter (18).Both the vector of mean values F and the covariance matrix C depend on the value of the parameter s.The case s = 0 trivially corresponds to F = F and C = C, while s > 0 leads to F ̸ = F and C ̸ = C.In order to quantify the deviation of F from F we introduce the quantity ∆ defined as The value of ∆ represents the average deviation of the new values F from the starting ones F measured with respect to the starting covariance.In other words, ∆ < 1 means that on average F deviates from F by less than one standard deviation.
As a further estimator of the deviation of F with respect to F , we introduce also the quantity η defined as The value of η can be smaller or larger than unity depending on whether | F i | is (on average) smaller or larger than In the same spirit, in order to estimate how the uncertainties of the new mean values, i.e.
ii , deviates (on average) from the starting ones σ i we introduce the quantity ϵ defined as In Table II we have collected the mean values F and the uncertainties σ corresponding to some values of the parameter s, which can be compared with the starting values F and σ.The values of the quantities ∆ (see Eq. ( 47)), η (see Eq. ( 48)), ϵ (see Eq. ( 49)) and of the percentage p of events passing the unitary filter (18) are also shown.
The values of p and ∆ increase for increasing s as expected, while both η and ϵ are found to be substantially constant.As the value of the parameter s varies from ≃ 0.5 to ≃ 7.5, the value of p ranges from ≃ 2% to ≃ 93% and the one of ∆ from ≃ 0.5 to ≃ 0.8, while η ≃ 1.02 and ϵ ≃ 0.7.43) and ( 45) for various values of the parameter s before the application of the unitary filter (18) (assuming Q 2 0 = 0).The second column show the (symmetrized) electroproduction data from Ref. [13] (see Table I).The values of the quantities ∆, η, ϵ (see Eqs. ( 47)-( 49)) and of the percentage p of events passing the unitary filter (18) are given in the last four rows.
Not only F ̸ = F and σ ̸ = σ, but also the correlation matrix of the data corresponding to the new PDF ( 46) is different from the starting one, namely This point is illustrated through the heat maps of Fig. 4 for various values of the parameter s before the application of the unitary filter (18).As expected, the correlations among first neighbors increases for s > 0. We observe a slight dependence of ρ ij on the value of the parameter s.
We now apply the unitary filter (18) and select only the subsets of events satisfying unitarity for each value of the parameter s.On such subsets we calculate the mean values, uncertainties and correlations for the form factor and the transverse susceptibility.In order to adopt a compact notation we will denote these quantities by F = {F i }, σ = {σ i } and ρ = {ρ ij } with ρ ij = C ij /(σ i σ j ) and i, j = 0, 1, ...(N + 1).
The mean values F and the uncertainties σ corresponding to the form factor points are shown in Table III for some values of the parameter s and compared with the starting values F and σ.The quantities ∆, η and ϵ are calculated using F and σ in Eqs. ( 47)- (49) and shown in Table III, as well as the size N sample = p • 10 5 of the subsets of events passing the unitary filter (18).It can be seen that the application of the unitary filter (18) leads to values {F and  II, but for the mean values F and uncertainties σ obtained using only the subsets of unitary events selected by the filter (18) for various values of the parameter s.The size N sample = p • 10 5 of these subsets is shown in the last row.
σ, which exhibit a weaker dependence on the parameter s with respect to the values F and σ obtained before the  (46) for various values of the parameter s before the application of the unitary filter (18).The case s = 0 corresponding to the initial correlation matrix ρij ≡ Cij/(σiσj) → δij is also shown.The labels 0, 1, ...7 correspond to the form factor points of Table I, while the label 8 corresponds to the non-perturbative result for the transverse susceptibility 4M 2 π χT (0).
application of the unitary filter.As the value of the parameter s varies from ≃ 0.5 to ≃ 7.5, the value of the quantity ∆ ranges only from ≃ 0.7 to ≃ 0.8, while both η ≃ 1.03 and ϵ ≃ 0.7 do not change significantly.Finally, also the correlation matrix ρ, obtained after the application of the unitary filter and shown in Fig. 5, changes only slightly with respect to the matrix ρ (see Fig. 4) obtained before the application of the unitary filter.
The values F of the form factor data are by construction normally distributed.After the application of the unitary filter (18) also the distributions of the new values F are substantially normal regardless the size of the sample of the unitary events, as shown in Fig. 6.Thus, we apply a simple iteration procedure to get rid off the residual dependence on the value of the parameter s.Namely, we apply the unitary filter (18) to the distributions of the values F and recalculate on the new subset of unitary events the mean values, uncertainties and correlations for the form factor and the transverse susceptibility.The percentage of events passing the unitary filter increases drastically and we iterate (few times) the previous steps until the percentage reaches ≃ 99%.This is done for various values of the parameter s and we find that the final set of values for F and for the covariance matrix C is almost totally independent on the starting value of s while keeping ∆ ≃ 0.8, η ≃ 1.03 and ϵ ≃ 0.7.We remind that the values obtained for ∆, η and ϵ means that on average the new central values F are larger w.r.t. the original ones F by few percent with changes lower than one standard deviation and that the new uncertainties σ are on average around 70% of the original ones  4, but for the correlation matrix ρ ij ≡ Cij/(σiσj) obtained after the application of the unitary filter (18) for various values of the parameter s (see Table III).
σ.We stress that F ̸ = F and C ̸ = C are direct consequences of the application of the unitary filter (18).
Using our final set of values for F and for the covariance matrix C we generate a sample of events for the input values of the form factor, which all satisfy the unitary filter (18).For the k-th event, corresponding to form factor values f ), we apply the DM formulae ( 23)-( 25) obtaining with with ϕ(z, Q  After summing over the sample we get the averages F L(U ) (Q 2 ), the standard deviations σ L(U ) (Q 2 ) and the correlation coefficient ρ LU (Q 2 ).Adopting a uniform distribution between F L (Q 2 ) and F U (Q 2 ) one finally obtains the following expressions for the em pion form factor The results obtained using Eqs.( 53)-( 54) starting from the electroproduction data of Ref. [13] are shown in Fig. 7, where they compared also with the CERN SPS data from Ref. [5].The latter ones are not used to construct the DM band.Nevertheless, the DM predictions at low Q 2 based on the electroproduction data at higher Q 2 deviates from the CERN data only within ∼ 1σ.
We close this Section by stressing that the unitary sampling procedure can be easily generalized to any set of  hadronic form factors, which must satisfy unitary bounds.

VI. UNITARY SAMPLING APPLIED TO THE CERN AND JLAB-π EXPERIMENTAL DATA
We now add the direct CERN data [5] to the electroproduction JLAB-π ones [13] obtaining a total of 52 data points distributed in the z-range [0.043, 0.70].We include also the data point F 0 = 1 at Q 2 0 = 0 (i.e.z 0 = 0) related to the charge conservation, obtaining a total of N + 1 = 53 data points.
The CERN data are uncorrelated with the electroproduction JLAB-π ones and, therefore, the covariance matrix is block diagonal, namely where, as in the previous Section, the JLAB-π covariance matrix C JLAB−π is diagonal, i.e. of the form given in Eq. ( 39), while for the CERN data it is necessary to include a normalization error δr = 0.45% [5] beyond the tabulated uncertainties σ i .We do that by using the following covariance matrix for the CERN data3 where F i ± σ i are the tabulated values of the pion form factor in Ref. [5].
We adopt a sample of 10 4 events generated according to the PDF (40) with the covariance matrix given by Eq. (55).As in Section IV we consider also a gaussian distribution for the non-perturbative transverse susceptibility 4M 2 π χ T (Q 2 0 = 0) = 0.00574 (10).This distribution is taken to be uncorrelated with those of the form factor points.Then, we calculate the values of χ DM (Q 2 0 = 0) (Eq.( 19)), which range from a minimum equal to ≈ 10 98 up to a maximum equal to ≈ 10 105 , i.e. extremely far from the 2-point bound 4M 2 π χ T (Q 2 0 = 0) = 5.74 (10) • 10 −3 .This is due to the fact that the kinematical coefficients d i (see Eq. ( 20)) have alternating signs with absolute values ranging from ∼ 10 14 to ∼ 5 • 10 55 .The huge cancellation occurring among the individual contributions to the r.h.s. of Eq. ( 19) can be handled using multiple arithmetical precision, achieved by adopting the software package MPFUN from Ref. [60] and using an adequate number of significant digits for the arithmetic operations.
The use of the PDF (46) with the modified mean values F and covariance matrix C, given respectively by Eqs. ( 45) and ( 43), allows to generate events fulfilling the unitary filter (18).As the parameter s varies from ≃ 20 to ≃ 45, the percentage p of unitary events ranges from ≃ 1% to ≃ 96%.Following the iterative procedure described in the previous Section we obtain a final set of values F for the form factor and C for the covariance matrix.We stress that both F and C are almost totally independent on the starting value of s, while the quantities ∆ ≃ 1.0, η ≃ 1.01 and ϵ ≃ 0.4 (see Eqs. ( 47)-( 49)) still remain at acceptable values.
The unitary form-factor values F turn out to be almost normally distributed.The correlation matrix corresponding to the final covariance matrix C is shown in Fig. 8 and compared with the initial correlation matrix of the input data C given by Eq. ( 55).The corresponding unitary band for the pion form factor, obtained using Eqs.( 53)-( 54 55)) and the final correlation matrix C, obtained using the unitary sampling (46) after 10 iterative steps.As in Fig. 4 the last label corresponds to the non-perturbative result for the transverse susceptibility 4M 2 π χT (0).
in Fig. 9 as the red band.It can be seen that the inclusion of the CERN data is very effective in producing a more precise band for the pion form factor at all Q 2 .Thanks to the DM method we can evaluate the slope of the pion form factor at Q 2 = 0 in a way completely independent of any (unitary) parameterization or explicit z-expansion.Using both the CERN and JLAB-π data (i.e. the red band of Fig. 9) we obtain for the pion charge radius ⟨r π ⟩ the result which is definitely larger than the PDG value (4) by ≈ 1.8σ.By neglecting the normalization error of the CERN data (i.e. by putting δr = 0 in Eq. ( 56)) we get a more precise result, ⟨r π ⟩ DM = 0.695 ± 0.014 fm, which differs from the PDG value (4) by ≈ 2.5σ.The DM result (57) differs also from the estimate ⟨r π ⟩ = 0.663±0.006fm made in Ref. [5] using the CERN data, but adopting a simple monopole Ansatz for the fitting function.We have checked that using the covariance matrix (56) a monopole fitting function leads to ⟨r π ⟩ = 0.656 ± 0.008 fm with a value of χ 2 /(d.o.f.) ≃ 1.0.However, when a dipole term is added to the monopole one, we get a quite different value of the pion charge radius, namely ⟨r π ⟩ = 0.699±0.024fm (again with a value of χ 2 /(d.o.f.) ≃ 1.0), which agrees much better with the DM result (57) for both the mean value and the uncertainty.These findings indicate clearly that the estimate of ⟨r π ⟩ made in Ref. [5] as well as those from Refs.[4,6] are plagued by a significative model dependence, so that they cannot be considered parameterization independent.
As already pointed out, the DM result (57) differ by ≈ 1.8σ from the determination of ⟨r π ⟩ obtained using the abundant and precise timelike e + e − data in Refs.[16,19], while exhibiting a much larger uncertainty.In order to clarify any possible significance of the above difference a significant improvement of the precision of the experimental data in the spacelike region is called for.
CERN data (sym.)JLAB-π data FIG.9.The red DM band (at 1σ level) corresponding to the input values F and C obtained after 10 iterative steps of the unitary sampling procedure applied to both the CERN SPS [5] (black circles) and the (symmetrized) electroproduction JLAB-π [13] (blue squares) experimental data.The blue DM band is the same as in Fig. 7 and corresponds to the use of only the electroproduction JLAB-π data.
In Fig. 10 our DM band is compared with the results of Ref. [16], based on a unitary analysis of both timelike e + e − and spacelike CERN data.A good overall agreement is observed up to Q 2 ≃ 1 (GeV/c) 2 (see lower panel), while a zoom in the low-Q 2 region (upper panel) shows that: • the use of the very precise and dense timelike e + e − data leads to the accurate result for the pion charge radius obtained in Ref. [16], namely ⟨r π ⟩ = 0.655 ± 0.003 fm; • the DM band is in better agreement with the spacelike CERN data w.r.t. to the results of the dispersive analysis of Ref. [16].In this respect, since the pion charge radius is correlated with the 2π contribution to the muon HVP term [61], it would be interesting to analyze the possible impact of the new e + e − → π + π − experimental data from the CMD-3 Collaboration [36] on such correlated quantities [62].
We close this Section by observing that the addition of the spacelike data of the F2 [4] and the (less precise) SELEX [6] experiments at FNAL do not change significantly the results shown in Fig. 9 and in Eq. (57).In particular, using as input the F2 [4], CERN [5] and JLAB-π [13] data sets (for a total of 66 data points including the absolute normalization at Q 2 = 0) we get which differs from the PDG value (4) by ≈ 1.6σ.Had we chosen the vector susceptibility 4M 2 π χ T (Q 2 0 = 0) to be equal to the value 0.00550 (4), which corresponds to the evaluation of Eq. ( 34) using the dispersive pion form factor |F V π (ω)| from Ref. [16] up to ω = 1 GeV, the value of the pion charge radius ⟨r π ⟩ DM would result to be 0.702 ± 0.026 fm, i.e. very close to Eq. ( 58).

VII. UNITARY BGL FIT
In this Section we perform a BGL analysis of the spacelike data after constructing a truncated z-expansion in which unitarity is built-in.
As known, in the BGL approach [26] the product of the pion form factor F V π (Q 2 ) times the kinematical function ϕ is analytic inside the unit circle |z| = 1 and, therefore, it can be expanded as CERN data (sym.)JLAB-π data FIG.10.The red DM band of Fig. 9 compared with the results of Ref. [16] (labelled CHS 2019), based on a unitary analysis of both timelike e + e − and spacelike CERN data [5] (black circles).In the lower panel also the (symmetrized) electroproduction JLAB-π [13] (blue squares) experimental data are shown.
where ϕ(z, Q ) is given by Eq. ( 27) and the coefficients a k (Q 2 0 ) (which are real and depend implicitly also on the choice of the auxiliary quantity t 0 appearing in the definition of the conformal variable ( 7)) are constrained by the unitary bound In order to lighten the notation, we will indicate hereafter the coefficients a k (Q 2 0 ) simply as a k .The z-expansion (59), truncated at some order N BGL , can be used as a fitting Ansatz to describe the spacelike data for the pion form factor, namely with the unitary bound given by NBGL k=0 The truncation introduces unavoidably a model dependence.In Ref. [26] it was proposed to look at the truncation error δF BGL π (Q 2 ), defined as which has un upper bound given by The weak point of the bound (64) on the truncation error (63) is that it may represent an upper limit on the difference between the true function ( 59) and the truncated fit ( 61) if and only if the coefficients a k with k = 0, 1, ...N BGL of the true function coincides exactly with those of the truncated fit.However, generally speaking, this is not guaranteed and, therefore, the truncation bound ( 64) is of limited use, particularly when the bound ( 62) is almost saturated.We now want to apply the truncated BGL fit (61) to the description of the spacelike pion data with the unitary bound (62) built-in.This can be achieved through a simple procedure based on a hyperspherical transformation (see, e.g., Ref. [63]) described in Appendix B. We generate a sample of 10 3 events F according to the PDF (40) using the direct CERN [5] and electroproduction JLAB-π [13] data with the covariance matrix C given by Eq. (55).At the same time, as in Sections IV and VI, we consider also a gaussian distribution for the non-perturbative transverse susceptibility 4M 2 π χ T (Q 2 0 = 0) = 0.00574 (10).This distribution is taken to be uncorrelated with those of the form factor points.
For each event we fit the input data using the BGL Ansatz (61) corresponding to a given truncation order N BGL with the unitarity bound (62) built-in through the hyperspherical procedure of Appendix B. Then, we minimize the reduced χ 2 r -variable given by obtaining the best unitary BGL fit for a given truncation order N BGL .We have considered values of N BGL between 2 and 10 and the corresponding results for the BGL parameters a k are shown in Table IV for N BGL = 2, 4, 6, 8, 10.Note that the value of the parameter a 0 is constrained by the absolute normalization F V π (Q 2 = 0) = 1 and by the value of the transverse susceptibility 4M 2 π χ T (Q 2 0 = 0) = 0.00574 (10), leading to a 0 = 0.190 ± 0.002 for any value of N BGL .The following comments are in order.
• As N BGL increases, the mean values and uncertainties of the coefficients a k with k ≲ 6 tend to remain stable.
• The coefficients a k for high-order monomials (k ≳ 6) have large uncertainties (up to ≈ 100%).
• For N BGL ≥ 4 the unitary bound ( 62) is almost saturated, while unitarity is strictly fulfilled (see in Table IV the maximum and average values of the parameter r 2 0 , defined in Eq. (B5)) .
• The values of the reduced χ 2 r -variable (65) are always consistent with unity.
• The distribution of the coefficients a k is approximately normal for low-order monomials (k ≲ 6), while significative deviations from a Gaussian distribution occur in the case of higher-order monomials (see Fig. 11) due mainly to the saturation of the unitary bound (62).Nevertheless, the form factors values are almost normally distributed, as shown in Fig. 12.
• For a given order N BGL of the truncation the parameters a k are generally anticorrelated (see Fig. 13).(61) for various values of the truncation order NBGL.The minimum (r 2 0 )min, maximum (r 2 0 )max and the average ⟨r 2 0 ⟩ values of the parameter r 2 0 representing the unitary bound (see Eq. (B5) for its definition) are shown together with the average value of the reduced χ 2 r -variable (65) corresponding to a sample of 10 3 events generated using the direct CERN [5] and electroproduction JLAB-π [13] data with the covariance matrix given by Eq. ( 55).
The bands for the pion form factors corresponding to the unitary BGL fits of Table IV are shown in Fig. 14 for various values of the truncation order N BGL .It can be seen that in the kinematical region covered by the CERN and JLAB-π data (i.e. for Q 2 ≲ 2.5 GeV 2 ) the results of the unitary BGL fit are stable against the order N BGL of the truncation, while at larger values of Q 2 the bands are largely unstable and no extrapolation is possible at least from N BGL ≤ 10.This finding clearly indicates the inadequacy of the estimate of the truncation error based on Eq. ( 64), since the saturation of the unitary bound (62) would imply a negligible truncation error that is not observed at all at large Q 2 .
Moreover, a closer look to Figs. 9 and 14 reveals that the BGL and DM bands differ by 1 ÷ 2σ in the kinematical region of the electroproduction JLAB-π data (i.e. for 0.35 ≤ Q 2 (GeV 2 ) ≤ 2.45).This observation will be discussed in a while.
The pion charge radius ⟨r π ⟩ corresponding to the BGL fit ( 61) is explicitly given by The results for ⟨r π ⟩ BGL obtained for the unitary BGL fits of Table IV are shown in Fig. 15 and exhibit a good convergence as a function of the truncation order N BGL .We get ⟨r π ⟩ BGL = 0.717 ± 0.044 fm, which is consistent with the DM result ( 57) with an uncertainty larger by a factor ≃ 1.5.Including also the spacelike data of the F2 experiment at FNAL [4] we obtain which differs from the PDG value (4) by ≃ 1.3σ.We mention that, if we adopt for the vector susceptibility 4M 2 π χ T (Q 2 0 = 0) the value 0.00550 (4), which corresponds to the evaluation of Eq. ( 34) using the dispersive pion form factor |F V π (ω)| from Ref. [16] up to ω = 1 GeV, the pion charge radius ⟨r π ⟩ BGL remains basically unchanged with respect to Eq. (67).
The instability of the unitary BGL fits of the pion form factor observed in Fig. 14 at Q 2 ≳ 2.5 GeV 2 and the larger uncertainty of the BGL result (67) w.r.t. the corresponding DM result (57) are connected to non-unitary effects present in the fitted CERN and JLAB-π experimental data.According to the DM method these data, corresponding to the input sets F and C for the pion form factor values and their covariance matrix, do not fulfill the unitary bound (11).The unitary sampling procedure, described in Section V, has allowed us to get in Section VI a new set of input data F and C fulfilling the unitary bound (11).IV.The black solid lines represent a gaussian fit of the various histograms.
Therefore, we apply the unitary BGL fit (61) with N BGL = 10 to the unitary set of input data F and C defined by the DM method, i.e. we replace in Eq. ( 65) the input data set F and C with the DM set F and C. We stress that only in this way unitarity is fulfilled both by the fitting function (by construction) and by the fitted data (by the unitary sampling procedure).The corresponding band for the pion form factor is shown in Fig. 16 as the green band.The comparison with the bands from Figs. 9 and 14 indicates very clearly that non-unitary effects are significantly present when the set F and C of input data are adopted regardless the fact that the BGL fit satisfies unitarity by construction.Such non-unitary effects produces not only a large instability of the BGL fits in the kinematical region not covered by the experimental data, but they can also affect significantly the fitting results in the electroproduction region.Finally, the pion charge radius corresponding to the green band of Fig. 16 is ⟨r π ⟩ = 0.707 ± 0.029 fm; both the central value and the uncertainty are now in nice agreement with the DM result (57).

VIII. UNITARY BCL APPROACH
An alternative z-expansion is the so-called BCL one, originally proposed in Ref. [27] to address the momentum dependence of the hadronic form factors describing the semileptonic B → πℓν ℓ decays.
In the case of the em pion form factor the BCL expansion is a direct z-expansion and its truncated form read as An interesting feature of the BCL approach is the inclusion of the analytic constraint at the annihilation threshold CERN data (sym.)JLAB-π data The pion charge radius ⟨rπ⟩ corresponding to the unitary BGL fit (61) for various values of the truncation order NBGL, applied to the CERN [5] and the (symmetrized) electroproduction JLAB-π [13] experimental data.The red band correspond to the DM result (57), while the black one to the PDG value (4) given in Ref. [2].z = −1.Indeed, angular momentum conservation requires that in the timelike region Im[F V π (ω)] ∝ (ω 2 − 4M 2 π ) 3/2 .In turn this implies that the real part of pion form factor should have a vanishing first derivative at the annihilation threshold, namely dF V π (z)/dz = 0 at z = −1.Such a constraint can be easily implemented in the truncated BCL  approach (at variance with the BGL one) by adding a further monomial term b N BCL +1 z N BCL +1 to Eq. ( 68), namely where the matrix B jk (Q ) and of the transverse susceptibility 4M 2 π χ T (Q 2 0 ) (see Ref. [27]).In Appendix C we extend the procedure used in the case of BGL fits to include the unitary condition (70) in a χ 2 -minimization fitting approach based on the BCL fit (69).We calculate also the matrix B jk (Q 2 0 ) for Q 2 0 = 0 and N BCL up to 10.We have applied the truncated BCL fit (69) to the description of the spacelike pion data with the unitary bound (70) built-in.As in the case of the unitary BGL fits performed in the previous Section, we have used a sample of 10 3 events F generated according to the PDF (40) using the direct CERN [5] and electroproduction JLAB-π [13] data with the covariance matrix C given by Eq. (55).A gaussian distribution for the non-perturbative transverse susceptibility 4M 2 π χ T (Q 2 0 = 0) = 0.00574 (10), uncorrelated with those of the form factor points, has been assumed.We have minimized the reduced χ 2 r -variable given by obtaining the best unitary BCL fit for a given truncation order N BCL .Note that the value of the parameter b 0 is constrained by the absolute normalization F V π (Q 2 = 0) = 1 to be equal to b 0 = 1 for any value of N BCL .We have considered values of the truncation order N BCL between 2 and 10, obtaining results similar to those of the BGL fits shown in Fig. 14 with a slightly better precision.In the kinematical region covered by the CERN and JLAB-π data (i.e. for Q 2 ≲ 2.5 GeV 2 ) the results of the unitary BCL fit are stable against the order of the truncation N BCL , while at larger values of Q 2 the bands are unstable and no extrapolation is possible at least from N BCL ≤ 10.We find that for N BCL ≥ 4 the unitary bound ( 70) is almost saturated.As in the case of the BGL fits performed in the previous Section, this finding implies the inadequacy of truncation errors based only on higher order terms in the BCL fit (69), since the saturation of the unitary bound would imply a negligible truncation error that is not observed at all at large Q 2 .
For the pion charge radius the results corresponding to the BCL fits exhibit a good convergence as a function of the truncation order N BCL , obtaining ⟨r π ⟩ BCL = 0.713 ± 0.031 fm, which agrees very well with the DM result (57).Including also the spacelike data of the F2 experiment at FNAL [4] we get which differs from the PDG value (4) by ≃ 1.8σ.As already observed in Sections VI and VII in the cases of the DM and BGL approaches, respectively, if we adopt for the vector susceptibility 4M 2 π χ T (Q 2 0 = 0) the value 0.00550 (4), which corresponds to the evaluation of Eq. ( 34) using the dispersive pion form factor |F V π (ω)| from Ref. [16] up to ω = 1 GeV, the pion charge radius ⟨r π ⟩ BCL remains basically unchanged with respect to Eq. (72).
In this Section we address the issue of the impact of the value of the auxiliary quantity Q 2 0 , at which the transverse susceptibility 4M 2 π χ T is evaluated, on the unitary filter (18) and on the corresponding DM band for the pion form factor.
The Q 2 0 -dependence of the l.h.s. of the inequality ( 18) is shown in Fig. 1 for increases and such a drop is mainly governed by the mass of the dominant ρ(775)-meson resonance in F V π (ω).Instead, the Q 2 0 -dependence of the r.h.s. of the inequality ( 18) is due to the kinematical function ϕ(z, Q 2 0 ), defined in Eq. ( 27) where z 0 is given by Eq. (22).By definition the function ϕ(z, Q ϕ(z, Q In this way, the DM band for the pion form factors can be obtained from Eqs. ( 23)- (25) simply by replacing χ T and ϕ with χ T and ϕ, namely and the unitary filter (18) becomes The Q 2 0 -dependencies of 4M ) obtained using for |F V π (ω)| the results of the dispersive analysis of the e + e − data available from Ref. [16] up to ω = 1 GeV and putting |F V π (ω)| = 0 for ω > 1 GeV.The dashed line corresponds to add to the data-driven pion form factor a power-law tail for ω > 1 GeV of the form which is very sensitive to the high-energy tail of |F V π (ω)|.
At low values of Q 2 (≲ 0.25 GeV 2 ) the sensitivity to the choice of Q 2 0 is quite limited and the more precise determination of the pion charge radius ⟨r π ⟩ is the one obtained at Q 2 0 = 0, presented in Section VI.In the Q 2 -region of few GeV 2 the DM band becomes more precise as Q 2 0 increases.In particular, the DM band at Q 2 0 = 0 may be consistent with negative values of the pion form factor for Q 2 ≈ 4 − 5 GeV 2 .This tendency is less pronounced as Q 2 0 increases and it is expected to disappear for Q 2 0 > 1 GeV 2 .Note that it is reasonable to exclude zeros in the pion form factor, particularly in the spacelike region, where in quantum mechanics the pion form factor is the Fourier transform of a charge distribution proportional to the square of the pion wave function 4 .
The main conclusion of this Section is that the choice of the value of Q 2 0 can have an impact on the DM predictions in the region of Q 2 ≈ few GeV 2 (and beyond), where the optimized choice is expected to be Q 2 0 ≈ few GeV 2 .On the contrary, for the pion charge radius ⟨r π ⟩, which represents an important quantity investigated in this work, the optimized choice is given by Q 2 0 = 0, as properly considered in Section VI.Thus, it would be valuable to obtain in the next future a reliable determination of the susceptibility 4M 2 π χ T (Q 2 0 ) for Q 2 0 ≳ 1 GeV 2 either coming from LQCD simulations or driven by analyses of timelike e + e − data.We close this Section by recalling that the precision of the unitary DM band for F V π (Q 2 ) depends also on the quantity and the Q 2 -range of the input data.In this respect, both the addition of experimental results for Q 2 up to ≈ 8.5 GeV 2 , planned at JLAB [65,66], or even up to ≈ 30 GeV 2 , expected at future facilities like Electron-Ion Colliders [67,68], and the inclusion of precise LQCD determinations of F V π (Q 2 ) at low and high Q 2 would be very valuable.

X. THE ONSET OF PERTURBATIVE QCD
The behavior of the pion form factor at large spacelike momentum transfer is predicted by the pQCD hard-scattering mechanism [32][33][34][35] to be where f π ≃ 130 MeV is the pion decay constant, is the running strong coupling at leading order (with N f being the number of active flavors and Λ QCD the QCD scale) and G(Q 2 ) describes the corrections due to the pre-asymptotic structure of the scale-dependent pion distribution amplitude, which is function of the light-front fraction of the pion's total momentum carried by a valence quark.The question at which energy scale the asymptotic behaviour (80) sets in has long been debated in literature and the answer is not trivial because of the presence of nonperturbative effects at intermediate values of Q 2 (see, e.g., the recent review in Ref. [69]).
However, the situation changes when Q 2 0 becomes sufficiently large.Indeed, from Eq. ( 74) for Q While the central value β(z) drops down as 1/Q 5/2 , i.e. faster than the pQCD prediction (80), the unitary width γ(z) does not.This means that at a certain value of Q 2 = Q 2 pQCD the pQCD prediction (80) may start to be within the unitary bounds β(z) ± γ(z) at 1σ level (for Q 2 0 >> Q 2 pQCD ).Such a value Q 2 pQCD provides an estimate of the energy scale at which the asymptotic behaviour (80) sets in and it is based only on unitarity and spacelike experimental data.
Note that in principle one should calculate the DM unitary bands for the pion form factor for increasing, but finite values of Q Since a precise estimate of the transverse susceptibility 4M 2 π χ T (Q 2 0 ) is not available for Q 2 0 > 1 GeV 2 , we limit ourselves to investigate the sensitivity of the DM unitary band to a range of possible values for the quantity 4M 2 π χ T (∞) given by Eq. (79).The first estimate is calculated using for |F V π (ω)| the results of Ref. [16] and cutting the integral in the r.h.s. of Eq. (79) at ω = 1 GeV.Due to the positivity of the integrand function in Eq. ( 79), such a value represents a lower bound to the transverse susceptibility, namely By considering for ω > 1 GeV a power-law tail of the form 4 an additional contribution equal to 0.009 ± 0.003 is obtained, leading to which corresponds to an increase of about 30%.Since the susceptibility 4M 2 π χ T (∞) is very sensitive to the high-energy tail of the em pion form factor, we consider for our sensitivity study a third higher value by applying to the result (83) a conservative factor equal to 2, obtaining 6 The DM bands for the quantity predicted by the DM method adopting the three estimates (82)-(84) for the transverse susceptibility 4M 2 π χ T (∞) in the limit Q 2 0 → ∞.The input data are the CERN [5] (black circles) and the (symmetrized) electroproduction JLAB-π [13] (blue squares) experimental data.The orange line represents the leading pQCD prediction (80) with G(Q 2 ) = 0 and ΛQCD ≃ 300 MeV [70].The black line corresponds to the monopole shape ) reproducing the PDG central value (4) for the pion charge radius.lower unitary bound (at 1σ level) turns out to be almost insensitive to the chosen value of 4M 2 π χ T (∞) and remains significantly much larger than the leading pQCD prediction (80) (with G(Q 2 ) = 0) at least up to Q 2 ∼ 100 GeV 2 .Therefore, the pre-asymptotic structure of the pion distribution amplitude is expected to produce significant effects on the pion form factor up to quite large values of Q 2 .This is in qualitative agreement with the findings of several estimates available in the literature based both on models [71,72] and on LQCD simulations [73][74][75].
We point out that the DM bands shown in Fig. 20, which we recall are based only on unitarity and spacelike experimental data, may provide important information on the scale dependence of the pion distribution amplitude. 6By assuming conservatively for ω > 1 GeV a tail of |F V π (ω)| proportional to the leading behaviour expected in pQCD for timelike momenta [35], i.e. ω −2 log −1 (ω 2 /Λ 2 QCD ) with Λ QCD ≃ 300 MeV [70], one gets an additional contribution to Eq. (82) equal to ≃ 0.04, yielding a value of 4M 2 π χ T (∞) not exceeding the one given in Eq. (84).
As noted in Section IX, the precision of the unitary DM band for F V π (Q 2 ) at large Q 2 can be improved by adding new experimental results, like those planned at JLAB [65,66] for Q 2 up to ≈ 8.5 GeV 2 .The projected precision of the forthcoming JLAB experimental data is shown in Fig. 21 by the green triangles and compared with the DM unitary band corresponding to the transverse susceptibility (84).
predicted by the DM method in the Q 2 -range up to 10 GeV 2 adopting the estimate (84) for the transverse susceptibility 4M 2 π χ T (∞) at Q 2 0 → ∞.The input data are the CERN [5] (black circles) and the (symmetrized) electroproduction JLAB-π [13] (blue squares) experimental data.The green triangles (fixed at an arbitrary value of 0.4) correspond to the projected precisions (statistical+systematic) expected in the experimental proposal [66] for the 12 GeV JLAB upgrade.The orange and black lines are the same as in Fig. 20.

XI. CONCLUSIONS
The experimental data on the em form factor of charged pions available at spacelike momenta have been analyzed using the DM approach [1], which describes the momentum dependence of hadronic form factors without introducing any explicit parameterization and includes properly the constraint coming from unitarity and analyticity.The latter one is given by a transverse vector susceptibility, which has been evaluated nonperturbatively from the results of lattice QCD simulations of suitable two-point correlation functions contributing to the Hadronic Vacuum Polarization term of the muon.
We have elucidated in detail the role played by the DM unitary filter (18), which allows to select the subset of input data that can be reproduced exactly by a unitary z-expansion.Since the unitary bound turns out to be extremely selective as the number of data points increases, we have develop a unitary sampling procedure, which allows to generate in an efficient way a distribution of values for the pion form factor satisfying unitarity for any value of the number of data points.Such a procedure can be generalized straightforwardly to any set of hadronic form factors, which must satisfy unitary bounds.
We have applied the unitary sampling method to the spacelike data from both the CERN SPS experiment [5] and the JLAB-π Collaboration [13] for a total of more than 50 data points.The pion charge radius has been determined in a completely model-independent way and consistently with the unitary bound.This is at variance with the results obtained in the experimental works [4][5][6], where the spacelike data have been fitted assuming a simple monopole Ansatz, which introduces a non-negligible model dependence.
The DM result is ⟨r π ⟩ DM = 0.703±0.027fm, which differs by ≃ 1.6 standard deviations from the latest PDG [2] value ⟨r π ⟩ P DG = 0.659 ± 0.004 fm, dominated by the very precise results of dispersive analyses of timelike data [16,19] coming from measurements of the cross section of the e + e − → π + π − process.In order to clarify any possible significance of such a difference it is crucial to improve significantly the precision of the experimental data in the spacelike region.We stress that the procedure described in this Appendix for constructing a unitary BGL fitting function can be applied to a generic hadronic form factor which should fulfill a unitary constraint.
Appendix C: Inclusion of the unitary constraint in truncated BCL fits In this Appendix we describe briefly the procedure that allows to span the space of the values of the BCL coefficients b k , appearing in the truncated BCL fit and satisfying the (truncated) unitary constraint Charge conservation requires that F V π (Q 2 = 0) = 1, which implies b 0 = 1 in Eq. (C1).Thus, the unitary constraint (C2) should be evaluated putting b 0 = 1.In the case the unitary sum exceeds unity, the coefficients b k with k = 1, 2, ...N BCL can be multiplied by a common factor chosen to ensure that the unitary sum is equal to unity, while b 0 is kept equal to unity.
Finally, as described in Section VIII, the inclusion of the analytical constraint at the annihilation threshold z = −1 corresponds to add in Eq. (C1) a further monomial term b NBCL+1 z NBCL+1 with the value of the coefficient b NBCL+1 fixed by those of the coefficients b k with k = 1, 2, ...N BCL .Such an addition requires to re-evaluate the unitary constraint (C2) with N BCL replaced by N BCL + 1.In the case the new unitary sum exceeds unity, the coefficients b k with k = 1, 2, ...N BCL + 1 can be multiplied by a common factor to ensure that the unitary sum is equal to unity, while the analytical constraint remains fulfilled.
The procedure described in this Appendix is not limited to the case of the em pion form factor, but it can be applied to a generic hadronic form factor which should fulfill a unitary constraint.

2 π χ T (Q 2 0
FIG.3.Histograms of the ratio χDM /0.00574 corresponding to the data susceptibility χDM (see Eq. (19)) obtained using samples of 10 5 events generated according to the multivariate distribution(46) for various values of the parameter s.The quantity 0.00574 represents the central value of the non-perturbative transverse susceptibility(38), obtained in Section III at Q

sFIG. 4 .
FIG.4.Heat maps representing the correlation matrix ρij ≡ Cij/( σi σj) corresponding to the unitary sampling(46) for various values of the parameter s before the application of the unitary filter(18).The case s = 0 corresponding to the initial correlation matrix ρij ≡ Cij/(σiσj) → δij is also shown.The labels 0, 1, ...7 correspond to the form factor points of TableI, while the label 8 corresponds to the non-perturbative result for the transverse susceptibility 4M 2 π χT (0).

sFIG. 5 .
FIG.5.The same as in Fig.4, but for the correlation matrix ρ ij ≡ Cij/(σiσj) obtained after the application of the unitary filter(18) for various values of the parameter s (see TableIII).

FFIG. 6 .
FIG.6.Distributions of the form factors values F obtained after the application of the unitary filter(18) for various values of the parameter s.The green, red, blue and orange histograms correspond respectively to Q 2 = 0.35, 0.70, 1.0, 2.45 GeV 2 .The black solid lines represent a gaussian fit of the various histograms.

FIG. 7 . 2 0
FIG.7.The DM band (at 1σ level) for the em pion form factor corresponding to the set of values F and the covariance matrix C obtained starting from the (symmetrized) electroproduction data analyzed by the JLAB-π Collaboration in Ref.[13], shown as the blue squares.The value Q 2 0 = 0 is assumed.The experimental data obtained at CERN SPS[5] (black circles) are shown just for comparison, but they are not used to construct the DM band (see text).

FIG. 8 .
FIG.8.Heat maps representing the initial correlation matrix of the CERN + JLAB-π data (see Eq. (55)) and the final correlation matrix C, obtained using the unitary sampling (46) after 10 iterative steps.As in Fig.4the last label corresponds to the non-perturbative result for the transverse susceptibility 4M 2 π χT (0).

10 FIG. 11 .
FIG. 11.Distributions of the parameter a2, a6, a8 and a10 of the unitary BGL fit(61) with NBGL = 10 shown in TableIV.The black solid lines represent a gaussian fit of the various histograms.

2 0
where b NBCL+1 = NBCL k=0 k b k (−) k−NBCL /(N BCL + 1).The unitary constraint for the coefficients b k is more involved w.r.t. the case of the BGL fit.It makes the coefficients b k dependent on Q and reads as NBCL+1 j,k=0 b

2 0
) does not know anything about meson resonances.Let us factorize out the term (1

2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 2 0 = 1
) is approximately flat at Q = 0, while it increases sizably at high values of z as Q increases.Thus, as Q increases, the quantity χ DM (Q ) increases, so that the r.h.s. of the DM filter (78) becomes more sensitive to the input data in the large Q 2 -region.Thus, from the above findings we can conclude that at Q = 0 the unitary filter (78) is dominated by the CERN data, which are more precise and dense w.r.t. to the JLAB-π data, while as Q increases the impact of the electroproduction data increases.As Q increases, both sides of the DM filter (78) increase.Whether this filter may lead to more precise form factor bands for Q > 0, can be established only by a direct numerical investigation.This has been done for three different values of Q not exceeding the limiting value Q GeV 2 , as discussed in Section III.The corresponding results are shown in Fig.19.

2 0 2 0 2 0 2 0
and then extrapolate such DM bands to the limit Q → ∞.In this way the analytic property of the kinematical function ϕ(z, Q ) is kept at each step of the calculation.In practice we have checked that at large Q 2 the DM band extrapolated to Q → ∞ can be obtained directly by considering ab initio the kinematical function ϕ(z, Q corresponding to the three choices (82)-(84) are shown in Fig. 20.The

FIG. 21 .
FIG. 21.The pion form factorQ 2 F V π (Q 2 )predicted by the DM method in the Q 2 -range up to 10 GeV 2 adopting the estimate (84) for the transverse susceptibility 4M 2 π χ T (∞) at Q

2 0 2 0 2 0 2 0 2 0 2 0
.[27], in the case of the em pion form factor, the matrix B jk (Q ) being the coefficients of the z-expansion of the kinematical function ϕ(z, Q ) divided by 4M 2 π χ T (Q the l.h.s. of Eq.(C4) represents an analytic, bounded function inside the unit disc |z| ≤ 1, the coefficients B |j−k| (Q ) can be evaluated numerically by truncating the sum over n in Eq. (C3) up to a finite order.We have performed such a calculation at Q = 0 using n ≤ 300 and 4M 2 π χ T (0) = 0.00574.The corresponding results for the first 12 coefficients B |j−k| (0) are B(0) = {+0.1470510,+0.0630217, −0.0396807, −0.0388320, −0.0148702, −0.0085370, −0.0056373, −0.0040269, −0.0030298, −0.0023662, −0.0019011, −0.0015618} , (C5) which allow to construct the matrix B jk (0) for N BCL ≤ 10.The distribution of values of the transverse susceptibility 4M 2 π χ T (0) can be taken exactly into account by dividing all the coefficients B |j−k| (0) given in Eq. (C5) by the common factor 4M 2 π χ T (0)/0.00574.Following the strategy described in the case of the unitary BGL fit in Appendix B and dropping for sake of simplicity the dependence upon Q , we introduce N BCL + 1 parameters r k with k = 0, 1, ...N BCL , whose values can vary by construction in the range [0, 1].Then, using the hyperspherical rotation defined by Eqs.(B3)-(B4), we transform the set of parameters r k into a set of (intermediate) coefficients a k satisfying the constraint b k can be obtained by observing that the matrix B jk is symmetric and positive definite, so that it can be diagonalized having only positive eigenvalues.Thus, we can obtain the N BCL + 1 coefficients b k from the set {a k } as where v m is the eigenvector of the matrix B corresponding to the eigenvalue λ m .It follows that NBCL j,k=0 b j B jk b k =

TABLE II .
Mean values F and uncertainties σ obtained from Eqs. (

TABLE III .
The same as in Table

TABLE IV .
Values of the parameters a k of the unitary BGL fit [13]unitary band (at 1σ level) obtained using the unitary BGL fit(61)for various values of the truncation order NBGL, applied to the CERN SPS[5](black circles) and the (symmetrized) electroproduction JLAB-π[13](blue squares) experimental data.
4(see Section III).The range of Q