Strong evidence of $\rho(1250)$ from a unitary multichannel reanalysis of elastic scattering data with crossing-symmetry constraints

An analysis is presented of elastic $P$-wave $\pi\pi$ phase shifts and inelasticities up to 2 GeV, aimed at identifying the corresponding $J^{PC}=1^{--}$ excited $\rho$ resonances and focusing on the $\rho(1250)$ vs. $\rho(1450)$ controversy. The approach employs an improved parametrization in terms of a manifestly unitary and analytic three-channel $S$-matrix with its complex-energy pole positions. The included channels are $\pi\pi$, $\rho2\pi$, and $\rho\rho$, the latter two being effective in the sense that they mimic several experimentally observed decay modes with nearby thresholds. In an alternative fit, the $\rho2\pi$ mode is replaced by $\omega\pi$, which is also an experimentally relevant channel. The improvement with respect to prior work amounts to the enforcement of maximum crossing symmetry through once-subtracted dispersion relations called GKPY equations. A separate analysis concerns the pion electromagnetic form factor, which again demonstrates the enormous importance of guaranteeing unitarity and analyticity when dealing with very broad and highly inelastic resonances. In the case of $\rho(1250)$ vs. $\rho(1450)$, the failure to do so is shown to give rise to an error in the predicted mass of about 170 MeV. A clear picture emerges from these analyses, identifying five vector $\rho$ states below 2 GeV, viz. $\rho(770)$, $\rho(1250)$, $\rho(1450)$, $\rho(1600)$, and $\rho(1800)$, with $\rho(1250)$ being indisputably the most important excited $\rho$ resonance. The stability of the fits as well as the imposition of unitarity, analyticity, and approximate crossing symmetry in the analyses lend very strong support to these assignments. The possibly far-reaching consequences for meson spectroscopy are discussed.


I. INTRODUCTION
The experimental status of meson resonances with masses ranging from 1 to 2 GeV is very poor. Many states expected from the quark model have not been observed so far, whereas several apparently normal resonances listed in the PDG tables [1] do not fit in with mainstream quark models like, for instance, the relativized meson model by Godfrey and Isgur (GI) [2]. In Ref. [3] some of the obvious discrepancies were briefly reviewed, e.g. concerning the many observed f 2 states (with J P C = 2 ++ ) to be contrasted with the much fewer ones predicted in the GI model. Another disagreement is the relatively low mass of the strange vector meson K (1410), the first radial excitation of K (892), which is predicted almost 200 MeV higher in the GI and similar quark models. On the other hand, ρ , the first radial excitation of ρ(770), is listed by the PDG as ρ(1450) [1], which is difficult to reconcile with a lighter K (1410), as the latter state contains one strange quark and one light quark instead of two light quarks. However, under the ρ(1450) entry in the PDG meson listings one finds a large variety of experimental observations, with a huge mass range of 1208-1624 MeV, also depending on the particular strong decay mode. As a matter of fact, there have been many indications of a lighter ρ , roughly in the range 1.25-1.3 GeV, which we shall examine in more detail in the next section.
The importance of more accurately knowing the ρ mass, and of course that of K , should not be underestimated. Quark models based on the usual Coulomb-plus-linear interquark potential, with a running strong coupling constant in the Coulombic part, predict increasing radial splittings for lighter quarks. Therefore, a ρ mass significantly lower than the value predicted in the GI and largely equivalent models would pose a serious challenge to such approaches. Now, it is true that the precise mass of a broad resonance like ρ depends on the way the corresponding scattering data are analyzed. In that respect, the usual Breit-Wigner (BW) parametrizations can be very unreliable, possibly leading to deviations of the order of 10 MeV for ρ(770) and even more than 100 MeV in the case of ρ [4]. The reason is that multi-BW parametrizations typically do not satisfy unitarity, which becomes a very serious issue in highly inelastic processes like those in which ρ is observed. Another problem, this time on the theoretical side, is the usual static approach to mesonic resonances in most quark models, treating them as manifestly stable bound states of a quark and an antiquark. Coupled-channel effects from meson loops, due to both open and closed meson-meson decay channels, may give rise to sizable mass shifts, alongside producing a large width. For instance, in the multichannel unitarized quark model of Ref. [5] a bare ρ mass of 1.48 GeV is lowered by about 190 MeV owing to coupled channels, resulting in a physical ρ mass of 1.29 GeV. However, in the latter model a completely different confinement interaction is employed, leaving doubts about the precise size of such unitarization effects in mainstream quark models. Moreover, also the ground states as well as the higher radial excitations may suffer con-siderable mass shifts [5], so that a refit of all parameters will have to be carried out when unitarizing any particular model.
A further important piece of information comes from recent lattice-QCD calculations. In Ref. [6] masses of light and strange hadrons were computed in an unquenched simulation, but with only (anti)quark interpolators included and so resulting in purely real spectra. Thus, the mass of the first radial excitation of K (892) was found to be slightly above 1.6 GeV. On the other hand, in Ref. [7] members of the same lattice collaboration study P -wave Kπ scattering in a simulation with both quark-antiquark and meson-meson interpolators. Here, they employ Lüscher's [8] method to extract Kπ phase shifts from discrete energy levels for different lattice sizes. The resulting phases are in good agreement with the data, including the K (892) mass and even its width in an extrapolation to the physical pion mass. The big surprise is the mass of the first radial excitation of K (892), coming out at (1.33 ± 0.02) GeV, that is, about 300 MeV lower than in the former lattice simulation, without meson-meson interpolators. Admittedly, the latter calculation amounts to an approximation, as only the Kπ decay channel is included, thus treating the K as a purely elastic resonance. Nevertheless, the importance of accounting for unitarity when doing spectroscopy in a quantitatively reliable fashion is unmistakable corroborated.
The present paper aims at clarifying the status of ρ and also the higher vector ρ excitations, by reanalyzing old data on ππ scattering, viz. elastic phase shifts and inelasticities up to about 2 GeV. The employed model of analysis is a manifestly unitary three-channel S-matrix parametrization, in which the complex pole positions of the different ρ resonances are explicitly included through generalized BW-type expressions. Moreover, (approximate) crossing symmetry is enforced by minimizing in the fits the difference between the experimental real parts of amplitudes and the theoretical ones resulting from dispersion relations. The three included channels are ππ, as well as the effective channels ρ2π and ρρ, with the latter ones mimicking 4π final states. For further details, see Sec. III.
The paper's organization is as follows. Section II extensively reviews the status of ρ(1250) vs. ρ(1450) in decades of literature. Section III describes the done Smatrix analyses of ππ phase shifts and inelasticities with the imposed crossing-symmetry constraints. In Sec. IV an analysis of the pion electromagnetic form factor further illustrates the necessity of a unitary and analytic approach to very broad inelastic resonances. Section V is devoted to a general discussion of the results and conclusions. The first time a ρ resonance was included in the PDG tables dates back to 1974 [9], with the entry in the datacard listings called ρ (1250). Its mass and width were listed as 1256 MeV and 130 MeV, respectively, from Ref. [10]. In the latter paper, evidence was found of two opposite-parity spin-1 ωπ resonances at about 1250 MeV inpp annihilations at rest. The vector state could correspond to the ρ and the pseudo-vector (J P C = 1 +− ) to what is nowadays called b 1 (1235) [1]. To our knowledge, though, the earliest indication of a possible ρ goes back to 1970 and was reported in Ref. [11], also cited in Ref. [9]. In this experiment, neutral bosons were observed in photoproduction on protons, including a resonance at about 1240 MeV with a width of roughly 100 MeV [11]. The authors tentatively identified this state with the pseudo-vector "B" (now b 1 (1235)) meson. However, not having determined the J P C quantum numbers, they concluded: "It is possible that this particle could be an as-yet-undiscovered vector meson" .
No further vector mesons were found in the energy region 1.3-2.0 GeV with a cross section larger than 5% of that of ρ 0 (770) (90% confidence level [11]).
In 1982 [20] and 1983 [5], 1982 two quark-model calculations [5,20] reported support for a ρ below 1.3 GeV, while also describing its small ππ cross section. Starting with Ref. [5], this above-mentioned unitarized multichannel quark model [5] was applied to vector and pseudoscalar mesons. The resulting P -wave ππ cross section reasonably reproduced the ρ(770) mass, width, and pole position, while also predicting a ρ pole with a real part of 1.29 GeV, though with a too small width. At this energy, the ππ cross section revealed a very small enhancement on top of the ρ(770) tail (see Ref. [5] , FIG. 4). Furthermore, the relativistic quark model of Ref. [20] predicted a ρ at about 1.22 GeV, explaining its small ππ width as a combined effect of the nodal structure of the radially excited ρ and the Lorentz-contracted wave functions of the outgoing two pions. Also of interest is the amplitude analysis of the coupled ωπ, ππ system in Ref. [21], which explained the difficulty to observe a ρ in ππ scattering between 1.1 and 1.3 GeV as being due to a small yet dominant inelastic background in that energy region.
The ρ(1250) entry in the PDG was maintained up till 1986 [22]. Things changed in 1987 with a combined anal-ysis [23] of two-pion and four-pion data from γp and e + e − processes, resulting in the postulation of two excited ρ resonances, with the masses 1.465 and 1.7 GeV. Curiously, no ρ below 1.3 GeV was even considered in the different fits to the data, in spite of its systematic inclusion in the PDG tables since 1974. As a matter of fact, no mention at all of such a resonance was made in Ref. [23]. On the other hand, the GI relativized quark model [2] mentioned above was explicitly cited, for having predicted the ρ at a mass of 1.45 GeV and the corresponding 1 3 D 1 state at 1.66 GeV, i.e., values very close to the ones found in these fits [23]. A key feature in the latter analysis is the authors' (unproven) conjecture, in the context of the vector-meson-dominance model for photoproduction, that the off-diagonal transitions 1 3 S 1 → 2 3 S 1 and 1 3 S 1 → 1 3 D 1 have cross sections comparable to the transitions 1 3 S 1 → 1 1 P 1 and 1 3 S 1 → 1 3 D 3 , which correspond to the diffractive photoproduction of the 1 +− "B(1250)" meson (now called b 1 (1235)) and 3 −− "g(1690)" (now ρ 3 (1690) [1]), respectively. The subsequent 1988 PDG edition [24] then included the new entries ρ(1450) and ρ(1700), while completely eliminating the ρ(1250) and oddly accommodating the ρ(1250) observation of Ref. [18] under ρ(1450). This state of affairs has remained unaltered so far [1].
One of the authors of Ref. [23] published several more papers on the ρ and related issues, which merit some further attention. The first of these appeared in 1991 [25], in a reaction to a new observation [26] of a ρ below 1.3 GeV, with mass 1266 MeV and width 166 MeV. These vector-resonance parameters resulted from a partial-wave amplitude analysis of the π + π − system observed with the LASS spectrometer at SLAC for the reaction K − p → π + π − Λ at 11 GeV. In Ref. [25] the authors claimed that "the interpretation of the LASS state at 1.27 GeV as a radial excitation of the ρ can in all probability be ruled out on the basis of its very small electromagnetic coupling" .
Note, however, that Ref. [26] did not make any statement about this coupling or the associated ρ(1270) → e + e − width. Instead, the authors of Ref. [25] carried out a BW fit including ρ(770), a "ρ x " with fixed m ρx = 1.27 GeV and Γ ρx = 0.17 GeV, a "ρ 1 " with fixed m ρ1 = 1.44 GeV and Γ ρ1 = 0.36 GeV, and a "ρ 2 " with mass and width to be adjusted to the data, resulting in m ρ2 = (1.73 ± 0.02) GeV and Γ ρ2 = (0.29 ± 0.07) GeV. The problem is that such a BW description of overlapping resonances violates S-matrix unitarity [27], which is all the more serious for extremely broad states, as in the present case. Moreover, the conclusion [25] that the found ρ x → e + e − width is too small was based on a comparison with GI model predictions [2], the very same ones that were claimed to be incompatible with a prior [23] analysis of ρ and ρ . Another curiosity in Ref. [25] is the following conclusion: "Thus there does not seem to be a strong case for the interpretation of the LASS state as an exotic" . This contrasts with several later collabora-tive works [28][29][30][31] of one of the authors of Ref. [25], which advocated the interpretation of the LASS "ρ x " resonance as a crypto-exotic four-quark state (having non-exotic quantum numbers). However, such an assignment would require [29] the existence of a narrow isoscalar partner state "ω x " at about 1.1 GeV, for which there is no experimental evidence. The very broad h 1 (1170) [1] is the only isoscalar J = 1 resonance with a nearby mass and it is easily interpreted as a normal unflavored 1 +− qq state.
Despite the dominant consensus on ρ at 1.45 GeV and the corresponding 1 3 D 1 state at 1.7 GeV, several later experiments and analyses contradict this picture. The OBELIX [32] spin-parity analysis ofpp → 2π + 2π − annihilations at very low momentum resulted in the clear identification of a ρ resonance with mass (1.282 ± 0.037) GeV and width (0.236 ± 0.036) GeV, i.e., values compatible with those [26] of the LASS collaboration. More recently, a combined 2-channel S-matrix and generalized BW analysis of P -wave ππ phase shifts and inelasticities up to 1.9 GeV was carried out [33], satisfying analyticity and multichannel unitarity. As a result, not only was a ρ(1250) firmly established for both methods of analysis, but evidence of higher ρ excitations was also found, viz. at roughly 1.6 and 1.9 GeV. Remarkably, a further state at about 1.45-1.47 GeV could be accommodated as well, though its inclusion in the fits turned out to be almost immaterial. A generalization of the mentioned 2-channel S-matrix parametrization to three channels in Ref. [34] largely confirmed these results, the most significant difference being the prediction of a higher ρ excitation at about 1.8 GeV instead of 1.9 GeV. Let us finish this discussion of the literature on ρ(1250) vs. ρ(1450) by quoting D. V. Bugg [35]: "It is not clear how to assign 1 − states to these [Regge] trajectories. The ρ(1450) does not seem to fit naturally as the first radial excitation of the ρ(770)." In the next section we shall outline in detail our present method of analysis, which amounts to a further improvement of the analyses performed in Refs. [33] and [34], by also imposing constraints from crossing symmetry.

III. ρ(1250) FROM AN ANALYSIS WITH CROSSING-SYMMETRY CONSTRAINTS
The most recent confirmation of ρ(1250), fully supported by physical and mathematical arguments, resulted from a dispersive analysis [36] of the amplitudes in three coupled P -wave decay channels, with a built-in crossing-symmetry condition in the ππ channel. A more detailed account of this work appeared in the PhD thesis of one of the present authors (VN) [37], from which we have selected the most important results and figures for the present paper and section.
Besides the ππ channel, for which experimental phase shifts and inelasticities were available, two additional, effective channels were included in the analysis, which should simulate the dominant three-and four-body decays of the ρ excitations. The problem here is that higher ρ resonances have many observed decay channels, which would be unfeasible to account for completely in our Smatrix approach, as it would lead to a proliferation of Riemann sheets and complex poles. Therefore, we are guided by the decay modes in the PDG listings [1], considerations from channel couplings, and the ππ phase shifts themselves. Thus, under the PDG ρ(1450) entry, we notice the modes a 1 (1260)π, h 1 (1170)π, π(1300)π, and ρ(ππ) S-wave , all included in the 4π decays. Now, a 1 (1260), h 1 (1170), and π(1300) all decay mostly to ρπ, so that we have three decay channels of an excited vector ρ state that do not lie very far apart and all lead to a quasi-final state of ρππ, with 4π being of course the true final state. And then there is the ρ(ππ) S-wave channel, which will naturally be dominated by the f 0 (500) (alias σ) resonance, with central threshold in the same energy ballpark. So we mimic these decays by including one ρ2π channel, with threshold at 1055 MeV. Of course, the mentioned four channels have central thresholds that lie 195-385 MeV [1] higher, but the involved resonances are extremely broad, so that opening our effective channel at the ρ2π threshold seems reasonable. As for the third effective channel to be included in our analysis, we note that an isovector vector state couples very strongly to a P -wave ρρ state, which is just a matter of recoupling coefficients of spin, isospin, and orbital angular momentum [5]. Indeed, under the PDG [1] ρ(1700) entry, we see the ρρ decay mode among the 'dominant' ρππ decays, besides the already considered a 1 (1260)π, h 1 (1170)π, and π(1300)π modes. Now, the central PDG threshold mass of the ρρ channel is 1550 MeV, but also the ρ is a relatively broad resonance, so that it is reasonable to open this channel at a somewhat lower energy. Thus, a value of 1512 MeV was obtained already in Ref. [34], upon fitting the P -wave ππ phase shifts, which is the threshold value we keep in the present analysis, too. Note that in Refs. [34,37] this effective channel was called ρσ, but we now prefer to call it ρρ. Clearly, the ρf 0 (500) channel will also contribute to the ρ(1700) decays just as to those of ρ(1450), in spite of this mode not being listed under the ρ(1700) PDG entry [1]. Nevertheless, the designation ρρ appears more appropriate for our third effective decay channel.
These three channels were included in the analyses in Refs. [34,36,37], the methods and results of which we discuss below. However, when considering ρ(1250) decays, special attention should be paid to the ωπ decay mode observed in several experiments [1]. So in order to have a more complete analysis, we now also do a fit replacing the ρ2π channel by ωπ, with threshold at 922 MeV. The corresponding results we will discuss below, right after presentation of methots used in our fits.
The simple phenomenological approach in this description is not based on some model with a phenomenological potential, but rather on a careful analysis of the S-matrix poles of certain resonances on different Riemann sheets of the complex energy plane. Among the free parameters in the amplitudes, fitted to the experimental data and to dispersion relations called GKPY equations [38], are the complex pole positions themselves, making the obtained results largely model-independent. This also allows to avoid problems from searching for poles by analytically continuing the amplitudes into the complex energy or momentum plane, as the pole postions are determined directly in the fits.
The amplitudes are fully unitary and analytical, viz. of the form where s is the effective two-pion mass squared and the S kl are S-matrix elements. For example, in the case of the ππ channel, such an element reads and expressions for other matrix elements are given in Eq. 6 of Ref. [34]. Phase shifts and inelasticities are denoted by δ 11 and η 11 . For simplicity, they will be just called δ and η further on in the text. The S-matrix factors S res and S bgr stand for resonant and background parts, respectively, while d res are the Jost functions, which contain all the dynamics of the interacting particles, both in individual channels and between them. The momenta in a given channel are denoted by k i and the uniformizing variable w is defined as where s 2 and s 3 are the thresholds of the second and third channel, respectively. The variable w transforms the eight-sheeted Riemann surface into a simpler complex plane. A resonance pole is given by with E r the resonance mass and Γ r its full width. So for s = s r , we have and the resonance contributions S res are defined as where M is the number of all poles. The background function has the form with Θ(s, s 2 ) the Heaviside function (= 1 for s > s 2 ) and where a and b are real numbers.
In the analysis of Ref. [34], (prior to Ref. [36]), the whole mathematical formalism described above was presented, and poles connected to a given resonance yet lying on different Riemann sheets were grouped into so-called clusters. These Riemann sheets depend on the number and type of analyzed channels. In the N -channel case, the S-matrix is defined on a 2 N -sheeted Riemann surface. Riemann sheets are numbered according to the signs of the imaginary parts of the relative momenta in all channels. So for three channels there are eight Riemann sheets on which the poles can lie, and they are num- The same amplitudes were then used in the later works [36,37]. However, three important changes were introduced: a) the restrictions on the movement of poles as a function of coupling between channels were removed; b) the tested amplitudes (in S, P , D, and F waves) were subjected to the limitation resulting from fulfilling crossing symmetry in the ππ channel up to about 1100 MeV; c) a threshold expansion with four parameters was used below about 640 MeV.
Removing the restrictions in point a) only lead to significant shifts of pole positions for two of them, associated with the ρ(770) cluster. They shift by several hundred MeV and thus produce very small phase shifts and inelasticity, typical for a weak background. It can be said that these become "background poles". Two other poles of ρ(770) and all the others of the higher ρ states shift by only a few MeV or less. The limitation in point b) is done by introducing in χ 2 total a component χ 2 CS corresponding to the mentioned crossing symmetry (CS). The imposition of this symmetry is controlled by where N = 26 is number of chosen energies (between the ππ threshold and 1100 MeV) at which χ 2 CS is calculated, and is fixed at 0.01 in order to make the size of χ 2 CS comparable to the other contributions to χ 2 total (i.e., the χ 2 data parts). Furthermore, ReA (in) (E) is the real part of the amplitude used to fit the data and the GKPY equations, while ReA (out) (E) is the same quantity yet calculated through the dispersion relations ds K II (s, s ) ImA I (in) (s ). (8) Here, C II st is the crossing matrix between ππ channels with isospin I and I , a I 0 is the S-wave scatteringlength vector for isospin I , and K II (s, s ) are the corresponding kernels for once-subtracted dispersion relations with imposed crossing symmetry. As demonstrated in Ref. [38], these dispersion relations produce significantly smaller errors in the computed amplitudes (actually in their real parts, from Eq. (8)) than the well-known Roy [39] dispersion relations with two subtractions. In practice, the integrals in Eq. (8) are done from the ππ threshold up to about 2 GeV, because data are lacking at higher energies and so-called driving terms are used thereabove. These driving terms have the same structure as the kernel, but are not related to the experimental input amplitudes A I (in) (s ). Their s and t dependence is given by Regge parametrization.
The amplitudes in Ref. [34] had a bad threshold behavior, i.e., they produced incorrect ππ scattering lengths. Nevertheless, this did not prevent obtaining very reasonable results for the resonance pole positions of the different ρs for lying far above the ππ threshold. However, when carrying out fits to the GKPY equations in Refs. [36,37], the threshold behavior of the amplitudes became very important, since the integrals in Eq. (8) start at threshold and the amplitudes there are the least suppressed by the dependence of the kernels on s (for explicit formulae of the kernels, see the Appendix of Ref. [38]).
In order to improve the near-threshold behavior of the amplitudes (point c) above) from Ref. [34], the original amplitude in Eqs. (1)-(6) is replaced by a polynomial below about 640 MeV (this value resulted from fits to the data and the GKPY equations). The polynomial is merely a generalized near-threshold expansion in powers of the pion momentum k, viz.
where a and b are just scattering length and effective range, respectively, which can be fixed or fitted to the data and to the dispersion relations. The parameters c and d are free in the fits to the data and to the GKPY equations, being used to smoothly match the phase shifts from the polynomial (i.e., their values and first derivatives) to the multichannel ones determined by Eqs. (1)-(6) at the matching energy of about 640 MeV. This value is still below the pole mass of ρ(770), so the effectiverange approximation can be used within these limits, as opposed to, for example, the S wave and the low-lying f 0 (500), in which case the matching energy must be well below 500 MeV. This way the ππ P wave is fitted simultaneously to the dispersion relations and the experimental data. The χ 2 total is defined as where χ 2 CS defined in Eq. (7) includes input from all 6 important partial waves (JI: S0, S2, P 1, D0, D2, and F 1) and χ 2 Data contains only phase shifts and inelasticities of the P 1 wave (hereafter just called P wave). Merely this wave is adjusted during the fits and changes to all ReA I(out) (s) in Eq. (8) are caused exclusively by modifications of this single wave. The free parameters are: mass and width of all P -wave resonances, i.e., for ρ(770) (8 parameters), ρ(1250) (16 parameters), ρ(1450) (8 parameters), ρ(1600) (8 parameters), and ρ(1800) (8 parameters), the matching energy, and the effective-range parameter a from the background part in Eq. (6) (b is fixed at −0.85 × 10 −4 to avoid violating unitarity in the inelasticities above 1.7 GeV). Thus, the maximum number of free parameters in the fits is 50. The parameters a and b in the polynomial defined in Eq. (9) are kept fixed at the values 0.0381 m −3 π and 0.00523 m −5 π , respectively. To avoid ending up in some local minima instead of the global one, the fits are performed sequentially with an increasing number of free parameters, viz. as follows: • 1st step: only matching energy, background parameter a, and ρ(770), i.e., 10 free parameters; • 2nd step: as in the 1st step, plus ρ(1250), i.e., 26 free parameters; • 3rd step: as in the 2nd step, plus ρ(1450), i.e., 34 free parameters; • 4th step: as in the 3rd step, plus ρ(1600), i.e., 42 free parameters; • 5th step: as in the 4th step, plus ρ(1800), i.e., 50 parameters.
Additionally, the fits are at each stage carried out repeatedly with a different number and values of the added resonance's initial parameters. The fitted parameters in each step serve as starting parameters in the next step. In the first step, also the matching energy and background parameter a are taken at different initial values. The total number of employed parameters (50) and their respective numbers for each resonance are the smallest ones that lead to good values of χ 2 total . Increasing them further does no longer give rise to a significant improvement in χ 2 total . As already mentioned above, a relevant decay mode not considered in our analysis so far is ωπ, which is included under both the ρ(1450) and ρ(1700) PDG [1] entries, with the corresponding resonance mass ranges of 1250-1624 MeV and 1550-1800 MeV, respectively. Our choice to do the main analysis with an effective ρ2π channel instead of the ωπ channel was based on the consideration that the effective one should be more important, as it is believed to account for several observed decay modes. Redoing our analysis while replacing the ρ2π channel by ωπ just amounts to changing the threshold value to 922 MeV, down from 1055 MeV, whereafter again fits to the experimental data and the GKPY equations are carried out. As a result, the masses of all resonances (i.e., the real parts of the poles) change only slightly: ρ(770) by +0.1 MeV, ρ(1250) by +9.7 MeV, ρ(1450) by +6.5 MeV, ρ(1600) by +2.4 MeV, and ρ(1800) by -7.5 MeV.
These small changes again show the stability and reliability of the obtained results. Moreover, comparing the quality of the two fits (values of χ 2 ) as presented in Table I, we see that they are essentially equivalent. The number of degrees of freedom (n.d.f.) in the fit equals 297, that is, 191 data points + 6×26 (26 energies for 6 partial waves) minus 50 free parameters. In the fit with the ωπ channel, the energy dependence of the phase shift and inelasticity does not undergo any significant qualitative modification and only small quantitative changes.  The following description focuses on our analysis with ρ2π as the second channel.
The introduction of the three modifications a), b), and c) does not cause as significant alterations in the P -wave as in the S-wave treated in Ref. [36], but it brings about several numerical changes. First of all, it is guaranteed that this new amplitude is not only unitary and analytic like in Ref. [34], but it also has the correct threshold behavior and fulfills the crossing-symmetry condition from the ππ threshold up to about 1100 MeV. This modified amplitude, when fitted to the data given in Ref. [36] as well as to the GKPY equations, results in a set of pole clusters from which these resonance poles, i.e., the ones having by far the largest effect on the whole amplitude, are presented in Table II.
In order to efficiently take into account the influence of all higher ρ decay channels not included in the d res (w) Jost function in Eq. (5), the simplest possible background is introduced (see Eq. (6)) and fitted to the data as well as the GKPY equations. As a result, a constant and small phase of almost -20 • and a smoothly increasing small inelasticity are obtained.
The results presented in Figs. 1-5 are based on the fits carried out in Refs. [36,37]. Figure 1 displays the ππ phase shifts and inelasticities fitted to the also shown experimental data. Inspecting this figure and Table I, we clearly see that the curves reproduce the data very well, all the way from the corresponding thresholds up  to almost 2 GeV, especially when considering the small errors of the phase shifts. Experimental data are from [40]- [45] for the phase shifts and from [40]- [43], [45] for inelasticity. Figure 2 shows the P -wave ππ phase shifts due to the individual resonances, corresponding to poles (all members of a cluster for a given resonance) on different Riemann sheets. Of course, only the full phase shift has the correct threshold behavior, given by a polynomial with fixed scattering length and effective range. As one would expect, ρ(770) has by far the largest influence on the overall phase, dominating the contributions of the ρ excitations. Moreover, the second most important resonance is clearly ρ(1250), whereas the smallest effect is due to ρ(1450). The visible yet rather insignificant kinks in the phase-shift and inelasticity curves right above 1.5 GeV correspond to the opening of the sharp effective ρρ threshold. Clearly, they do not affect the quality of the fits at all. Nevertheless, from a theoretical point of view it would be desirable to somehow smear out this threshold so as to account for the ρ width, and the same for the ρ2π channel. An empirical way to do this was formulated in Ref. [46], by allowing for complex masses in the final state. The resulting violation of Smatrix unitarity was then corrected by redefining S with the help of a factorization valid for an arbitrary complex symmetric matrix. This may be the subject of future work along the lines of the present analysis.
To properly assess the contribution of individual resonances to the full amplitude, it is very clarifying to compute it before and after removing those resonances.  [45] for the phase shifts and from [40]- [43], [45] for inelasticity.
nates between 1.0 GeV and 1.5 GeV, being comparable to that of ρ(1600) and ρ(1800) thereabove. The ρ(1450) contribution is quite small over the entire tested energy range. Figure 4 shows the ππ inelasticity η for the full amplitude and also the individual resonances. One can see very well that even below 1.5 GeV (near the ρρ threshold) inelasticity due to the ρ(1250) amplitude significantly differs from 1 and together with part from that of ρ(770) almost completely determines the inelasticity of the full amplitude. Contributions from ρ(1450) and ρ(1600) largely cancel each other between the ρ2π and ρρ thresholds. Above roughly 1.5 GeV, ρ(1800) determines the energy dependence of η almost entirely, interfering with the still large but already rather unstructured ρ(1250) part comparable in size to that of ρ(770). The contribution of ρ(1450) to η is very small above 1.5 GeV. As expected, it only has a minor maximum at about 1.4 GeV. The significant drop in the full inelasticity at about 1.6 GeV is mostly determined by ρ(1600), after the opening of the ρρ channel. The role of the background in building η is small, showing a slow and smooth rise. Just as in the case of the phase shift, the energy dependence of the inelasticity for the amplitude without a given resonance, i.e., by omitting all poles associated with it on the different Riemann sheets, is very informative. In Fig. 5, we see that removing ρ(1250) would cause the largest change (after that caused by ρ(770)) to the inelasticity curve as compared to the one due to the full amplitude. Similarly, a significant modification would be caused by leaving out ρ(1600) or ρ(1800), but only around 1600 MeV or thereabove, respectively. Fi- nally, also here we observe that without ρ(1450) there would only be a modest change to η, over a relatively small energy region below 1.5 GeV, having little effect on the shape of the inelasticity curve. Considering the relationship between resonances and the inelasticity, one should refer to Eqs. (2)-(5). According to unitarity, the Jost functions d res in Eq. (5) are constructed in such a way that their ratio in the unitary S-matrix (Eqs. (1,2)) has a modulus equal to one in the whole elastic region. This is due to the full symmetry between the poles and the zeros there. In the inelastic region (s > s 2 ), this symmetry is automatically removed (w is no longer purely imaginary) and the moduli of numerators and denominators are not equal anymore. The large number of complex poles (characteristic of such a mul- tiresonance analysis in a three-channel approach) needed to describe the data and to meet the crossing-symmetry condition naturally leads to various interferences among all these poles, resulting in a total inelasticity consistent with unitarity and the conditions imposed on the fits. Of course, the moduli of the d * res (−w * )/d res (w) ratios corresponding to individual resonances (i.e., clusters of corresponding poles) do not have to fulfill the unitarity condition of being smaller than 1, as can be seen in Fig. 4. The same is true for the tiny contribution to the inelasticity due to the background term defined in Eq. (6). Only mutual interferences among all resonances (poles) and the background produce the physical (unitary) result. The sometimes reported unitary inelasticity of a particular resonance corresponds to such a physical result (full η in Fig. 4), but is limited and calculated in the energy range selected for a given resonance.

IV. ρ(1250) FROM AN ANALYSIS OF THE PION ELECTROMAGNETIC FORM FACTOR
Vector isovector mesons below 2 GeV play a very important role in, for example, the determination of the pion electromagnetic (EM) form factor. When analyzing cross sections of e + e − → π + π − production, this form factor F EM,I=1 π (s) appears explicitly and contains information on the dynamics of all these mesons: where the pion "velocity" β π (s) = s−4m 2 π s , R is the amplitude for ρ-ω mixing interference, and the phase φ = arctan ω . An analysis in Ref. [4] compared two different approaches to determining F (EM)I=1 π (s). The first one was based on the popular Gounaris-Sakurai (GS) [47] model constructed by assuming that, for a wide energy range of the elastic region up to 1 GeV, the P -wave isovector ππ scattering phase shift satisfies a two-parameter effectiverange formula of the Chew-Mandelstam type, i.e., where q is the pion momentum in the CM system and h(s) is a simple logarithmic function of q and s. The pion EM form factor is then given by Using the fact that at δ = π/2 the left-hand side of Eq. (12) vanishes and, comparing with a Breit-Wigner distribution formula, the first derivative of the phase shift can be given by 1/m ρ Γ ρ , one can express F (EM)I=1 π (s) directly in terms of the ρ(770) mass and width.
The other approach employed in Ref. [4] was based on a very simple unitary and analytic (UA) formula with two symmetric poles and zeroes representing a resonance. The left-hand cut was simulated with one pole (and a symmetric zero). In this way, it was shown [4] that fits to e + e − → π + π − experimental cross sections and to elastic ππ phase shifts (from GKPY equation) made independently using the GS model and the UA approach give very different results for the ρ(770) resonance parameters, especially its mass. As one can see in Table III, the mass difference is much larger than the estimated errors, which are similar in size to those in the PDG tables [1].  More significant differences were found in Ref. [4] for the two higher ρ states, in fits to cross sections up to pion momentum squared s = 9 GeV 2 . The authors of this analysis pointed out that a generalization of the GS model above the inelastic threshold done in many experimental works "is without any deeper physical background, as the original G.-S. model for the ρ 0 meson contribution was constructed from the Pwave isoscalar ππ scattering phase shift given by the generalized effective-range formula of the Chew-Mandelstam type, which is evidently valid only in the elastic region." Nevertheless, to check the quantitative differences between the GS model and the UA approach above 1 GeV, the authors did as in some experimental analyses, i.e., they carried out fits to the data with three ρ states, viz. ρ(770), ρ(1450), and ρ(1700) using the GS model. Then they compared the results to those obtained from their enhanced unitary and analytic amplitude with poles as degrees of freedom. The results are presented in Table IV and compared to those from the PDG [1]. It is worth noting that the number of free parameters in the UA analysis is smaller (11) than in the GS model (14). And despite the fact that the values of χ 2 per degree of freedom are 1.84 (UA) and 0.98 (GS), the results of the UA analysis are much more realistic. In particular, we should draw attention to the 171 MeV lower ρ(1450) mass found in the UA approach as compared to the GS model, which difference is much larger than the reported errors in Table IV. Similarly, significant differences can be found in the literature for the mass of the very same resonance, which can often be easily explained by comparing their values found in e.g. BW and simple "pole" approaches. In a BW amplitude, the mass of a resonance is determined by the energy M BW at which the phase shift passes 90 • (M BW = 2 k 2 BW + m 2 for two interacting equal-mass particles). The unitary "pole" amplitude (like the UA one used above and in Ref [4]) has, for one single resonance, two symmetric poles (and corresponding zeros) at k r = ±a − ib. Then, the phase shift is given by δ = arctan 2bk a 2 +b 2 −k 2 and it is clear that the value δ = 90 • is attained for k = k BW (i.e., 2 √ a 2 + m 2 = M BW ). This difference gets larger according as b increases. For ρ(770) difference between M BW and "pole" mass M r defined by real part of 2 k 2 r + m 2 is about 8 MeV, which explains the discrepancy seen in Table III and found in many other theoretical and experimental analyses [1]. For ρ(1450) this discrepancy is several dozen MeV, which is less than the differences seen in Table IV. However, above 1 GeV one is dealing with a few very broad and highly inelastic ρ resonances and a simple reasoning like for ρ(770) is completely insufficient, requiring to also account for other phenomena typical of resonance interference. It is worth noting here that by making fits using the BW and "pole" approaches and comparing M BW and M r for the ρ(770) one can additionally request that the decay width be the same in both approaches. Then the difference between these masses is about 5 MeV.
A very good example of such an analysis, which is both qualitative and quantitative in explaining the phase-shift behavior around 1250 MeV and the determination of the pole position in the amplitude, is Ref. [21]. In this work entitled "Why is the ρ (1250) not Observed in the ππ Scattering ?", an M -matrix parametrization of the ππ and ωπ partial amplitudes is analyzed and compared with the parametrization as a sum of the inelastic resonance term and the background amplitude. A most important qualitative conclusion is that even a small inelastic background around 1300 MeV can completely hide a ρ(1250) in the ππ channel, leading to a non-resonant behavior there of the ππ phase shifts. The quantitative results of this analysis confirm this conclusion: in order to describe the experimental data very well, a background of only a few degrees around 1300 MeV is sufficient, the real part of the pole then comes out at about 1220 MeV, the width is roughly 320 MeV, and Γ ππ /Γ ωπ ≈ 0.15, the latter ratio being in reasonable agreement with the available experimental data [1].   [4] above the inelastic threshold up to pion momentum squared s = 9 GeV 2 . For a better comparison, the values from the PDG tables [1] are also given.

V. CONCLUSIONS
The results of our combined analyses unmistakable demonstrate the necessity to include a ρ resonance at about 1.26 GeV. The stability of the fitted pole positions as well as the manifest fulfillment of multichannel unitarity and optimized crossing symmetry in our approach lend strong support to the reliability of our excited ρ states, including the ones at about 1.42 GeV, 1.60 GeV, and 1.78 GeV. Straightforward spectroscopic arguments then impose the following quark-model assignments: ρ(1250)/2 3 S 1 , ρ(1450)/1 3 D 1 , ρ(1600)/3 3 S 1 , and ρ(1800)/2 3 D 1 . Confirmation of these four states, which were already found in a previous analysis [34], poses serious problems to mainstream quark models, unless at least ρ(1250) is interpreted as a crypto-exotic tetraquark state, for which there is no experimental or theoretical support (also see the discussion of "ω x " above). A ρ at 1.25 GeV is very hard to reconcile with the GI [2] model and similar ones, based on a Coulomb-plus-linear confining potential with a running strong coupling constant α s (q 2 ). The only way out would be to consider the GI ρ(1450) a "bare" quark-antiquark bound state that upon unitarization turns into the physical ρ(1250) resonance, similarly as in the unitarized quark model of Ref. [5]. However, the latter model also predicts mass shifts for all other states, especially the ground state ρ(770), besides employing a completely different confinement mechanism. Therefore, the complete spectrum including unitarization effects would have to be computed again in the GI and similar models, after refitting the parameters. But even more seriously, it is practically inconceivable that the 3 3 S 1 state in the GI model at 2.0 GeV could be lowered by unitarization to 1.6 GeV. The problem is that the radial splitting between the first and second excitation for mesons with light quarks in the GI model is larger than 500 MeV [2]. This gives rise to huge discrepancies with the observed spectra [1] not only for vector states, but also for tensors like the f 2 [3,34], as already mentioned in Sec. I. Therefore, one must either present very convincing arguments why some of the resonances identified in the present and previous [33,34] analyses should be interpeted as crypto-exotic states, or consider the possibility that the Coulomb-plus-linear confining potential with a running coupling constant is inadequate, at least in the way it is usually implemented in quark models.
Regardless of these considerations, we believe to have made a convincing case for ρ(1250), which should finally be rehabilitated in the PDG tables with a separate entry. Further expertimental analyses that do respect multichannel S-matrix unitarity would be most welcome, of course, besides realistic model and lattice calculations accounting for unitarization effects.
The problem with new experimental data is that they will most likely result from production processes and not elastic scattering, making their direct inclusion in an analysis as described in the present paper problematic. A possible way out would be applying the formalism for relating production and scattering amplitudes as outlined in Ref. [48] to our multichannel S-matrix approach. Also this will be a topic of future research.