Global parameterization of $\pi \pi$ scattering up to 2 ${\mathrm {\,GeV}}$

Abstract

We provide global parameterizations of $\pi \pi \rightarrow \pi \pi$ scattering S0 and P partial waves up to roughly 2 GeV for phenomenological use. These parameterizations describe the output and uncertainties of previous partial-wave dispersive analyses of $\pi \pi \rightarrow \pi \pi$, both in the real axis up to 1.12 ${\mathrm {\,GeV}}$ and in the complex plane within their applicability region, while also fulfilling forward dispersion relations up to 1.43 ${\mathrm {\,GeV}}$. Above that energy we just describe the available experimental data. Moreover, the analytic continuations of these global parameterizations also describe accurately the dispersive determinations of the $\sigma /f_0(500)$, $f_0(980)$ and $\rho (770)$ pole parameters.

1 Introduction

The unprecedented high statistics on hadronic observables attained at experiments like LHCb, Belle or BaBar require rigorous and precise parameterizations of final state interactions. Future Hadronic facilities (FAIR, PANDA, etc.) will be even more demanding. One of the most needed parameterizations is that of $\pi \pi \rightarrow \pi \pi$ scattering, since two or more pions appear very frequently as final products of many hadronic interactions. In addition, a renewed interest on $\pi \pi \rightarrow \pi \pi$ scattering is coming from lattice calculations, which have been recently able to obtain scattering partial waves with almost realistic masses [1].

Data on $\pi \pi \rightarrow \pi \pi$ scattering were obtained in the 70s [2,3,4,5,6] indirectly from the $\pi N\rightarrow \pi \pi N'$ reaction. Unfortunately, this technique gave rise to several conflicting data sets. Thus, for decades, crude models were enough to describe such data. The exception is the very low-energy region, both in the experimental and theoretical fronts. On the one hand, there are very precise data below the kaon mass coming from $K_{l4}$ decays [7, 8], particularly after the NA48/2 results [9]. On the other hand, Chiral Perturbation Theory (ChPT) [10, 11] provides a systematic and accurate low-energy expansion in terms of pion masses and momenta.

However, for most phenomenological applications the low-energy region is not enough, since the production of pions is generically more copious around resonances. ChPT can be successfully extended to the resonance region by means of dispersion relations [12,13,14,15], usually called Unitarized ChPT. Different versions or approximations of this method generate or reconstruct all resonances in $\pi \pi \rightarrow \pi \pi$ up to 1.2 GeV: the $\sigma /f_0(500)$ the $\rho (770)$ and the $f_0(980)$ [16,17,18,19,20] and even in $\pi K$ scattering. However, the prize to pay is the loss of a controlled systematic expansion, which hinders the calculation of uncertainties, and the length of the analytic expressions once one deals with coupled channels above $K \bar{K}$ threshold. Above 1.2 GeV one can introduce by hand other resonances, yielding a successful description of data [21], although with the same caveats as before and with expressions even more elaborated. Nevertheless, the interest of these unitarized approaches is that they can connect with QCD through the chiral parameters. Moreover, they provide a good semi-quantitative approximation, including values of resonance poles, which are much better than the usual description of two-pion interactions in terms of simple popular models, like the superposition of simple resonant shapes, Breit-Wigner formulas in different versions, isobar models, etc...For a recent review on modern $\pi \pi$ scattering dispersive determinations, UChPT and other models for the $\sigma /f_0(500)$ see [22].

The interest of those naive popular models is, on the one hand, their simplicity, since for most applications just the phase and the elasticity functions are needed, not an elaborated model of the interactions with other channels. On the other hand, they can be fairly reasonable for narrow isolated resonances, like the $\rho (770)$. However, such simple models provide an incorrect description of the scalar-isoscalar partial wave. In particular this is the case of the very broad $\sigma /f_0(500)$ pole and its interplay with the very narrow $f_0(980)$, together with the singularity structure in terms of cuts in the complex s plane. Actually, the rescattering of two pions in this channel is frequently described with some sort of Breit-Wigner parameterization for the $\sigma /f_0(500)$, which might be able to describe a wide bump in the data, but fails to describe the chiral constraints in the threshold region as well as the phase shift in the whole $\sigma /f_0(500)$ region. Recall that by Watson’s Theorem [23] any strong elastic rescattering of two pions must have the very same phase of the $\pi \pi \rightarrow \pi \pi$ partial-wave with the same isospin and angular momentum.

In general, modern Hadron Physics demands more precise and model-independent meson-meson scattering parameterizations. This has been achieved over the last two decades by means of dispersion relations aiming at precision, not only for $\pi \pi$ [24,25,26,27,28,29,30,31,32], but also for $\pi N$ [33, 34] or $\pi K$ scattering [35]. Unfortunately, we have found that, for the hadron community, these dispersive results, either obtained numerically from complicated integral equations or parameterized by piecewise functions, are not always so easy to implement or do not cover a sufficiently large energy region. Hence, the purpose of this work is to provide relatively simple and ready-to-use parameterizations of the phase and elasticity of the scalar-isoscalar and vector $\pi \pi \rightarrow \pi \pi$ scattering partial waves up to almost 2 GeV. They will be consistent with data globally from threshold up to approximately 2 GeV, and with the dispersive analysis in [28], which extends up to 1.43 GeV in the real axis. Moreover, we will impose that these parameterizations will provide a simple analytic continuation to the complex plane, consistent with the dispersive representation and the values for the pole positions and residues of the $\sigma /f_0(500)$, $\rho (770)$ and $f_0(980)$ resonances found in [36]. In addition, both the dispersive results for the threshold and subthreshold regions are also described by the global parameterization, consistently with the scattering lengths, slope parameters and S0 wave Adler zero values obtained in [28].

2 The input to be described

As we already commented, there are several $\pi \pi \rightarrow \pi \pi$ scattering data sets extending up to almost 2 GeV [2,3,4,5,6]. These are customarily given in terms of partial waves $t^I_\ell$ of definite isospin I and angular momentum $\ell$. We will also use the spectroscopic notation where the $\ell =0,1,2,3\ldots$ waves are referred to as S, P, D, F... waves, followed by their isospin. Unfortunately, all those data sets are often incompatible from one another and, moreover, simple fits to each separated set or to averaged data sets do not satisfy well dispersion relations [28, 37, 38, 38,39,40]. Nevertheless, it is possible to use dispersion relations as constraints to obtain a Constrained Fit to Data (CFD) [28] that still describes the $\pi \pi \rightarrow \pi \pi$ data on partial-waves but satisfies dispersion relations within uncertainties. Furthermore, the CFD fulfills the normality requirements of the residual distribution [41], hence ensuring that the standard approach for error propagation can be used. This CFD parameterization will thus be part of our input.

One might wonder why not using directly this CFD parameterization and why in this work we are trying to obtain another one. After all, this parameterization has become quite popular and it has been used in many phenomenological applications. There are several reasons.

First, the dispersion relations used in [28] are of two kinds and they were applied up to different energies, always below 2 GeV. One kind consists of a set of Forward Dispersion Relations, which were studied up to 1.43 GeV. These equations are rather simple, but unfortunately cannot be extended to the complex plane in search for poles. They are only useful as constraints on the real axis. The other kind consists of two sets of partial-wave dispersion relations, usually referred to as Roy equations [24, 25, 30, 31, 42,43,44,45,46,47,48] (with two subtractions) and GKPY equations [28] (with one subtraction). The former are more stringent in the low-energy region and the latter in the resonance region. Unfortunately, these partial-wave equations are limited to 1.12 GeV, although they can be rigorously continued to the complex plane in search for resonance poles. The existence of these different energy regions motivated the authors in [28] to describe the data with a piecewise parameterization, which in principle cannot be extended rigorously to the complex plane. Therefore, our first aim is to provide a rather simple but global analytic parameterization, with realistic uncertainties, that can be used from $s=0$ to 1.43 GeV. Thus, it will mimic the CFD piecewise parameterization in the real axis above the $\pi \pi$ threshold, which will be used as the first of our inputs to be described.

Second, the $\sigma /f_0(500)$ pole lies so deep in the complex plane that a careful dispersive determination is needed in order to extract its precise parameters rigorously [30, 36, 49]. Using the CFD parameterization as input in the GKPY equations, it was obtained numerically that its pole lies at $\sqrt{s_\sigma }=(457^{+14}_{-13})-i(279^{+11}_{-7}){\mathrm {\,MeV}}$ with a residue $\vert g\vert =3.57^{+0.11}_{-0.13}$. Now, the low-energy piece of the CFD parameterization [28] was constructed as a conformal expansion valid up to 850 MeV, which lies within the elastic $\pi \pi \rightarrow \pi \pi$ region. This CFD conformal piece can be continued to the complex plane finding $\sqrt{s_\sigma }=(474\pm 6)-i(254\pm 4){\mathrm {\,MeV}}$, which is fairly close, but it is not the pole obtained from the dispersive representation. This discrepancy does not improve when one includes further constraints in the real axis. Namely, even if the CFD conformal parameterization is extended up to the $K\bar{K}$ threshold to take into account the $f_0(980)$ effect or to the subthreshold region, in order to describe the dispersive value for the Adler zero, one still finds sizable discrepancies with the GKPY pole result. This problem was observed time ago [50,51,52,53]; arbitrarily small changes in the real axis input data may lead to indefinitely large variations for the analytic continuation to the complex plane. This illustrates how trying to obtain the $\sigma /f_0(500)$ pole from a data fit that only reaches 850 MeV is not precise enough. Actually, the effects of the $f_0(980)$ and other singularities, like the left hand cut, are significant at this level of precision. Hence, our second aim is to provide a simple analytic parameterization that reproduces simultaneously the dispersive poles of the $\sigma /f_0(500)$ and $f_0(980)$ and their interference. Thus, the numerical results of the GKPY dispersion relations in the complex plane, including the numerical values of the $\sigma /f_0(500)$ and $f_0(980)$ poles, and in the subthreshold region will be the second input to be described. For the P-wave we will proceed similarly, but just for the $\rho (770)$ pole.

Finally, the CFD parameterization and the dispersive data analysis from which it was obtained only reach 1.43 GeV, but there are more data up to almost 2 GeV. However, the data at those high energies have many well-known caveats. Some of them were already discussed in detail in [22, 54] and in appendix C of [37], but we summarize them here. First, in that energy region we have to rely on a single scattering experiment, the CERN-Munich Collaboration, so that systematic uncertainties relative to other experiments are not available. Nevertheless there is some information on partial waves for the final $\pi ^0\pi ^0$ system above 0.8 GeV from the GAMS experiment on $\pi ^-p\rightarrow \pi ^0\pi ^0 n$ [55], which will be of interest later. Second, the CERN-Munich collaboration has many different solutions for the $\pi \pi$ scattering partial waves. Of these, the most popular one for the S0-wave is the one published in 1973 [2], also called “solution b” in the collaboration compilation of Grayer et al. [3]. This solution is also consistent with a later reanalysis with polarized targets [5]. In addition, there is the “solution (− − −)”, which was the most favored in the 1975 collaboration reanalysis [4] and the most used solution for the P-wave. Note that both “b” and “(− − −)” solutions are compatible with one another below 1.43 GeV. Other solutions for both waves were already disfavored in that very same analysis. Third, both solutions have caveats. On the one hand, the inelastic contribution to all hadronic cross sections are expected to dominate over the elastic ones (something that has been verified for $\pi N$, KN and NN scattering). However, this is not the case of “solution b”. It is hard to understand why this should be different for pions. On the other hand, if the inelasticity is large, then it can be proved theoretically [56, 57] that the solution in terms of phase and elasticity is not unique. “Solution b” is an example of an almost elastic case and “solution (− − −)” of a strong inelastic effect. Finally, the very same convergence of the partial-wave expansion could be questioned at those energies, since around 1.7 GeV the F-wave is as large as the P-wave, the D0 wave as large as the S0 and the D2 actually larger than the S2. Furthermore, the solution (− + −) has been recently revisited and slightly modified by one of the authors of the CERN-Munich Collaboration [58]. He argues that it is still consistent if the I=2 wave is not considered elastic, finding some qualitative agreement with the GAMS experiment. In particular they both show hints of the $f_0(1500)$ resonance, in contrast to the other solutions. Note that this (− + −) solution is also compatible with the “b” and (− − −) solutions below 1.43 GeV.

Therefore, in view of the caveats above, we have extended our fits beyond 1.43 GeV using as our third source of input different sets of data. Namely: (I) the data of [2, 3, 5], which reaches up to 1.9 ${\mathrm {\,GeV}}$, to obtain a “solution I”, or II) the (− − −) data of [4], which reaches up to 1.8 ${\mathrm {\,GeV}}$, to obtain a “solution II”, or III) the updated solution (− + −) in [58] to obtain a “solution III”. Below 1.43 GeV the input is the same for all solutions and they agree within uncertainties. As a technical remark, we have ensured that the central value and the first derivative of both the phase and elasticity are continuous at the matching point, which is chosen at 1.4 GeV to avoid fitting the very end of the CFD parameterization. In any case, one should keep in mind that none of these solutions has been checked against dispersion relations above 1.43 GeV. Thus, beyond that energy they should be considered purely phenomenological data fits.

3 Analytic parameterizations

In this section we present the parameterizations used for the scalar-isoscalar and vector $\pi \pi \rightarrow \pi \pi$ partial waves. Let us first note that below the $K \bar{K}$ threshold the process will be considered elastic and hence it will be uniquely characterized by its phase shift $\delta ^I_\ell (s)$, as it is customary, through the following definition:

$$\begin{aligned} t^I_\ell (s)=\frac{\hat{t}^I_ \ell (s)}{\sigma (s)}=\frac{e^{i\delta ^I_\ell (s)}\sin \delta ^I_\ell (s)}{\sigma (s)}= \frac{1}{\sigma (s)}\frac{1}{\cot \delta ^I_\ell (s)-i}, \end{aligned}$$

(1)

where $\sigma _i(s)= 2 q_i(s)/\sqrt{s}=\sqrt{1-4m_i^2/s}$ is the two-body phase space, $q_i(s)$ being the CM momentum of a particle with mass $m_i$. For brevity we will suppress the subindex in the pion case and write $\sigma (s)\equiv \sigma _\pi (s)$. The elastic region will be described with conformal maps for both the S and P waves.

Let us also recall the standard inelastic partial-wave representation

$$\begin{aligned} t^I_\ell (s)=\frac{\eta ^I_\ell (s) \mathrm {e}^{2i \delta ^I_\ell (s)}-1}{2i \sigma (s)}, \end{aligned}$$

(2)

where the elasticity parameter $\eta ^I_\ell (s)$ and phase shift $\delta ^I_\ell (s)$ will be described by two independent real functions.

Let us now present separately the parameterizations we have used to describe the two partial waves of interest for this work.

3.1 S0-wave parameterization

As explained above, our parameterizations will be consistent with the dispersive data analysis of [28], which extends up to 1.43 ${\mathrm {\,GeV}}$. Above that we will only provide three phenomenological fits to three sets of incompatible data, carefully matched to our parameterizations below. Let us discuss both regions separately.

3.1.1 S0-wave parameterization below 1.4 ${\mathrm {\,GeV}}$

The $\sigma /f_0(500)$ and $f_0(980)$ resonances dominate the behavior of the S0 partial wave in this region. The somewhat controversial $f_0(1370)$ couples very weakly to two pions and its effect in this region can be treated as background. For our purposes it is important to remark that the $\sigma /f_0(500)$ has an associated pole very deep in the complex plane that produces a wide structure increasing monotonously from threshold up to roughly 900 MeV, reaching a phase-shift of $90^{\circ }$ around 800 MeV, as seen in Fig. 1. It is known [40] that the $\sigma /f_0(500)$ pole can be generated in the S0 partial wave by a simple truncated conformal expansion, that we will call $t^0_{0,\text {conf}}(s)$. However, above 900 MeV the $f_0(980)$ pole adds a further sharp increase that makes the phase larger than $200^{\circ }$ right below the $K\bar{K}$ threshold. The phase then keeps growing slower but monotonously up to 2 GeV.

It is worth noticing that the interplay between the $\sigma /f_0(500)$ and $f_0(980)$ poles produces a sharp dip in the modulus of the amplitude and the elasticity right above the $K\bar{K}$ threshold. In order to describe the $f_0(980)$ effects accurately and consistently with the dispersive results, we will factorize in the S matrix the $f_0(980)$ shape separately from the conformal expansion that contains the $\sigma /f_0(500)$ pole. In other words, we will use $S^0_0=S^0_{f_0}S^0_{0,\text {conf}}$, where

$$\begin{aligned} S^0_0= & {} 1+2i\sigma (s)t^0_0\,,\nonumber \\ S^0_{0,\text {conf}}= & {} 1+2i\sigma (s)t^0_{0,\text {conf}}\,,\nonumber \\ S^0_{f_0}= & {} 1+2i\sigma (s)t^0_{f_0}\,. \end{aligned}$$

(3)

This factorization ensures elastic unitarity for the S0 wave, i.e., $\vert S^0_0\vert =1$, when both the conformal and the $f_0(980)$ contributions fulfill elastic unitarity independently, i.e., $\vert S^0_{0,\text {conf}}\vert =\vert S^0_{f0}\vert =1$. This will be the case in the elastic region below the $\bar{K}K$ threshold, $s<4m_K^2$.

For our purposes, we are interested in the amplitude partial-wave

$$\begin{aligned} t^0_0(s)=t^0_{0,\text {conf}}(s)+t^0_{f_0}(s) +2i\sigma (s)t^0_{0,\text {conf}}\,(s) t^0_{f_0}(s). \end{aligned}$$

(4)

Now, the conformal factor of the partial wave is constructed by analogy to the elastic formulation in Eq. (1)

$$\begin{aligned} t^0_{0,\text {conf}}(s)= \frac{1}{\sigma (s)}\frac{1}{\Phi ^0_0(s)-i}, \quad s<1.4 {\mathrm {\,GeV}} \end{aligned}$$

(5)

where, building on [40]

$$\begin{aligned} \Phi ^0_0(s)=\frac{\sqrt{s}}{2q(s)}\frac{m_\pi ^2}{s-z_0^2/2}\left( \frac{z_0^2}{m_\pi \sqrt{s}}+\sum _{n=0}^N{B_n \omega (s)^n}\right) . \end{aligned}$$

(6)

Let us remark that $\Phi ^0_0(s)=\cot \delta ^0_0(s)$ in the elastic region $s\le 4m_K^2$, implying $\vert S^0_{0,\text {conf}}\vert =1$. The $s-z_0^2/2$ denominator provides the so-called Adler zero at $s_{Adler}=z_0^2/2$ required by chiral symmetry [59]. For $z_0=m_\pi$ one would recover the Current-Algebra result, namely the leading order ChPT value. However, for us it will be a free parameter, fundamental to describe the subthreshold region. As we will see, it comes out from the fits consistent with the dispersive evaluation, which in turn is consistent with higher order ChPT evaluations (see [28, 60]). Note also that if we did not include the $\sim 1/\sqrt{s}$ term added to the conformal series then, for the values of $B_n$ obtained from the fit, we would find in the unphysical region that $\Phi ^0_0(s_{Adler})\rightarrow +i\infty$ and $\Phi ^0_0(0)\rightarrow i0^+$, so that necessarily $\Phi _0^0(s)=i$ somewhere between those two subthreshold points. This would yield a spurious pole, i.e., a ghost, in the partial wave. As explained in [40] these ghosts are mostly harmless and have little relevance in the fit quality and the pole positions, but as a matter of principle it is better to remove them. Actually, as explained there, adding the $1/\sqrt{s}$ term to the conformal series simply shifts the value $\Phi ^0_0(0)$ so that the spurious pole disappears. The $1/\sqrt{s}$ cut does not introduce any additional singularity in the partial wave since its square-root cut falls right on the left cut. At the same time this term is suppressed in the physical region and thus barely affects the fit and resonance pole positions. As shown below, for this wave it will be enough to set $N=5$ to obtain a good overall $\chi ^2/d.o.f.$ in the elastic region.

The conformal variable is defined as

$$\begin{aligned} \omega (s)=\frac{\sqrt{s}-\alpha \sqrt{s_0-s}}{\sqrt{s}+\alpha \sqrt{s_0-s}}, \end{aligned}$$

(7)

where $s_0$ corresponds to the highest value of s where the expansion is real and then $\alpha$ sets the center of the conformal expansion. We have found in practice that the S0 wave is more conveniently described if the conformal expansion, by becoming imaginary, introduces some inelasticity above the $K \bar{K}$ threshold [50]). Thus, we choose $s_0= 4 m_K^2$ with $\alpha =1$ for simplicity, so that the expansion center lies near 0.7 ${\mathrm {\,GeV}}$. Hence, between $K\bar{K}$ threshold and 1.4 ${\mathrm {\,GeV}}$, the $\Phi ^0_0(s)$ function will be complex, which effectively introduces an inelasticity. We tried otherwise, with a higher $s_0$, so that no inelasticity would come from the conformal factor. We then found that the $t^0_{f_0}(s)$ would require many extra parameters with strong correlations. In addition these parameters would have huge and unnatural scale differences among themselves.

Actually, with our choice of conformal expansion we can use a simple and intuitive functional form for $t^0_{f_0}$, inspired by the expression used in [50]. Namely

$$\begin{aligned} t^0_{f_0}(s)=\frac{s \,G }{M-s-\bar{J} (s,m_\pi )\,s\, G -\bar{J}(s,m_K) m_K^2 f(s)}, \end{aligned}$$

(8)

which ensures elastic unitarity for $s<4m_K^2$. The $\bar{J}$ two-meson loop functions, or Chew-Mandelstam functions [61], provide the unitarity cut above each threshold and are defined as

$$\begin{aligned} \bar{J}(s,m_i)=\frac{2}{\pi }+\frac{\sigma _i(s)}{\pi }\log \left( \frac{\sigma _i(s)-1}{\sigma _i(s)+1}\right) . \end{aligned}$$

(9)

Note that the constant G in the numerator in Eq. (8) is multiplied by s in order to cancel the phase-space pole at $s=0$ in Eq. (4). In addition, this factor suppresses the inelastic contribution at low energies, hence ensuring that the low-energy region is dominated by the conformal parameterization. In principle, f(s) could be any real analytic function and for convenience we will build it as an expansion of Chebyshev polynomials. The main advantage of this procedure is the low correlation among parameters, which will provide a more realistic description of the uncertainties. Note that the expansion variable will not be s, but a linear transformation that maps the $[2 m_K,1.5{\mathrm {\,GeV}} ]$ energy region into the $[-1,1]$ segment, where Chebyshev polynomials are orthogonal. This variable is:

$$\begin{aligned} \omega _1(s)=2\frac{\sqrt{s}-2m_K}{1.5{\mathrm {\,GeV}}-2m_K}-1. \end{aligned}$$

(10)

Thus, the real function f(s) will be expanded as

$$\begin{aligned} f(s)=\sum _{i=0}^N K_i x_i(\omega _1(s)), \end{aligned}$$

(11)

where $x_i$ is the Chebyshev polynomial of order i and $K_i$ are fitting constants. In practice it is enough to set $N=3$ and so we will do. This function also suppresses the $f_0(980)$ contribution far from its nominal mass.

As a matter of fact, one can get an acceptable $\chi ^2/d.o.f.$ using Eq. (8) to fit the dispersive results in the real axis and the complex plane around the $K\bar{K}$ threshold. However, when so doing the $f_0(980)$ pole position does not come out at the precise dispersive value given in [28]. For this reason we will impose the dispersive value of its pole position in the fit, by fixing the G and M constants.

Let us then briefly recall how to reach the second Riemann sheet in search for poles. According to the S-matrix unitary relation $S S^{\dagger }=\mathbb {1}$ and taking into account the Schwartz reflection symmetry, $t^I_\ell (s+i \epsilon )=t^I_\ell (s-i \epsilon )^*$, then the partial wave in the second Riemann sheet $t_{l}^{(2),I}$ is algebraically related to itself in the first Riemann sheet by

$$\begin{aligned} t^{(2), I}_ \ell (s)=\frac{t^{I}_\ell (s)}{1+2i \sigma ^{}(s) t^{I}_\ell (s)}, \end{aligned}$$

(12)

where the $\sigma (s)$ determination is chosen so that $\sigma (s^*)=-\sigma (s)^*$ to ensure the Schwartz reflection symmetry

$$\begin{aligned} t^{(2), I}_ \ell (s^*)=t^{(2), I}_ \ell (s)^*. \end{aligned}$$

(13)

As a consequence, a pole $s_p=s_R+i s_I$ in the second Riemann sheet implies that a zero of the S matrix exists also in the first Riemann sheet at $s_p$. This imposes two constraints on Eq. (8), which allow us to fix G and M as follows:

$$\begin{aligned} M&=\frac{(f_I J_{K R}+f_R J_{K I}) (s_I (J_{\pi I}-2 \sigma _R)-s_R (J_{\pi R}+2 \sigma _I))}{d}\nonumber \\&\quad +\frac{(J_{\pi I}-2 \sigma _R) \left( s_I^2+s_R^2\right) }{d}-(f_I J_{K I}-f_R J_{K R}), \\ G&=-\frac{f_I J_{K R}+f_R J_{K I}+s_I}{d},\nonumber \end{aligned}$$

(14)

where we have defined the constants

$$\begin{aligned} f(s_p)&=f_R+i\,f_I, \nonumber \\ \bar{J}(s_p,m_K)&=J_{K R}+i\,J_{K I}, \nonumber \\ \bar{J}(s_p,m_\pi )&=J_{\pi R}+i\,J_{\pi I}, \nonumber \\ \sigma (s_p)&=\sigma _R+i\,\sigma _I, \nonumber \\ d&=J_{\pi I} s_R+J_{\pi R} s_I+2 (\sigma _I s_I- s_R \sigma _R), \end{aligned}$$

(15)

and $s_{p}=s_R+i s_I$ corresponds to the $f_0(980)$ pole position, which is therefore a parameter to be varied within its uncertainties in our formulas.

In summary, for the scalar wave below 1.4${\mathrm {\,GeV}}$ we will use Eq. (4) with $t^0_{f_0}(s)$ defined in Eqs. (8), (9) and  (14), whereas $t^0_{0,\text {conf}}(s)$, containing the $\sigma /f_0(500)$ pole, is defined in Eqs. (5) and (6), which above the $K \bar{K}$ threshold gives and additional contribution to the inelasticity besides that of $t^0_{f_0}$.

3.1.2 S0-wave parameterization above 1.4 ${\mathrm {\,GeV}}$

As we have emphasized repeatedly, from 1.43 ${\mathrm {\,GeV}}$ there are no dispersive data analyses and, besides, the data can be grouped into three inconsistent data sets. However, we are frequently asked if we could extend our parameterization beyond 1.43 ${\mathrm {\,GeV}}$. Thus, we will provide simple phenomenological fits to the three sets of data that we will match to our formulas below 1.4 GeV so that the whole parameterization and its derivative are continuous. For this we need the values at $s_m=(1.4 {\mathrm {\,GeV}})^2$ of the phase shift, the elasticity and their derivatives with respect to the energy squared, denoted with a prime. These inputs will be taken from the parameterizations below 1.4 ${\mathrm {\,GeV}}$.

To reduce the number of parameters, we will make use again of Chebyshev polynomials to describe the phase shift above $s_m$, namely

$$\begin{aligned} \delta ^0_0(s)&=\delta ^0_0(s_m)+\Delta ^0_0 [x_1(\omega _2(s))+1]+d_0 [x_2(\omega _2(s))-1]. \nonumber \\&\quad +d_1 [x_3(\omega _2(s))+1]+d_2 [x_4(\omega _2(s))-1],\nonumber \\ \Delta ^0_0&=\frac{\delta ^{\prime \,0}_0(s_m)}{\omega '_2(s_m)}-d_0 x'_2(-1) -d_1 x'_3(-1)-d_2 x'_4(-1), \nonumber \\&=\frac{\delta ^{\prime \,0}_0(s_m)}{\omega '_2(s_m)}+4d_0-9d_1+16d_2. \end{aligned}$$

(16)

The variable for the Chebyshev polynomials now is:

$$\begin{aligned} \omega _2(s)&=2\frac{\sqrt{s}-\sqrt{s_m}}{2 {\mathrm {\,GeV}}-\sqrt{s_m}}-1. \end{aligned}$$

The presence of $\delta ^0_0(s_m)$ in Eq. (16) ensures the continuity of the parameterization and the value of $\Delta ^0_0$ the continuity of the derivative. As we will see when fitting the data, we will need just two Chebyshev polynomials to get a good $\chi ^2/d.o.f.$ for solutions I and II, leaving just one free parameter, $d_0$, for those fits. However, three free parameters are needed to obtain an acceptable $\chi ^2/d.o.f.$ for solution III. In all cases the value of $\delta '^0_0$ will be kept fixed to the central value when calculating uncertainties. This means that although the derivative is continuous, the uncertainties on the derivative at that point might have a small kink. Otherwise the error band becomes unrealistically large.

Concerning the elasticity function, it will be fitted through an exponential function with a negative exponent, to ensure $0\le \eta ^0_0 \le 1$. We have found that in the case of the S-wave, Chebyshev polynomials in this exponent produce unwanted oscillations. Thus we will use a simple phenomenological expansion in powers of $Q(s)\equiv q(s)/q_m-1$, where $q_m=q(s_m)$. In practice we have found that five terms are needed at most to obtain a good fit. Explicitly:

$$\begin{aligned} \eta ^0_0(s)=\exp \left[ -\left( \sum _{k=0}^{4} \epsilon _k Q(s)^k\right) ^2\right] . \end{aligned}$$

(17)

Continuity at the matching point fixes

$$\begin{aligned} \epsilon _0=\sqrt{-\log (\eta ^0_0(s_m))}, \end{aligned}$$

(18)

and then the continuity of the derivative fixes

$$\begin{aligned} \epsilon _1=-\frac{4q_m^2}{\epsilon _0} \frac{{\eta ^0_0}'(s_m)}{{\eta ^0_0}(s_m)}. \end{aligned}$$

(19)

In practice only three free parameters will be needed at most to obtain an acceptable $\chi ^2/d.o.f.$, and just one will be enough for solution I. Note that the logarithm in Eq. (18) appears in the constants needed for the smooth matching, but it does not introduce any spurious analytic structure. As with the phase, now $\eta ^{0\,\prime }_0$ will be kept fixed to its central value when calculating uncertainties.

In summary, the S0-wave high-energy parameterization has six free parameters at most, but we will see that when fitting data solution II needs only four and solution I just two, setting to zero the remaining ones. Recall that up to 1.4 GeV all three solutions are compatible among themselves.

3.2 P-wave parameterization

The $\pi \pi$-scattering P-wave is completely dominated by the $\rho (770)$ meson, which is customarily described using simple resonance models, like variations of Breit-Wigner parameterizations. In many cases, this is fair enough. However, even though the $\rho (770)$ is usually considered as the prototype of narrow resonance, its width is relatively large compared to its mass, which explains that the $\rho$-meson shape cannot be fully described with precision using a simple Breit-Wigner function or within an Isobar Model, but requires additional shape parameters [62, 63]. Let us also recall that the $\rho (770)$ is the main player of vector meson dominance. Actually, it saturates the most common hadronic observables, like, for instance, the hadronic total cross section $\sigma (e^+e^-\rightarrow \text {hadrons})$, which implies applications well beyond low-energy meson physics. Thus, given its relevance for Hadron Physics, we will provide in this section an analytic parameterization to describe the $\pi \pi$ vector-isovector channel up to approximately 2 GeV.

This wave is much simpler than the S0, since the inelasticity sets in at much higher energies and is much smaller than for the S0 wave. Actually, there is no need to factorize explicitly any resonance pole in the inelastic region as we did for the $f_0(980)$, but for convenience we will use a similar analytic formalism to introduce the small inelasticity above $K\bar{K}$ threshold. This will allow us to continue analytically our partial wave to the complex plane and mimic the dispersive results of [28] within the Lehmann ellipse. Once more, we will separate the energy regions below and above 1.4 ${\mathrm {\,GeV}}$, because the latter is not tested against the dispersive representation and has inconsistent data sets that will be fitted separately later.

3.2.1 P-wave parameterization below 1.4 ${\mathrm {\,GeV}}$

Thus, as we did for the S0 wave, we build our partial wave as $S^1_1=S^1_{1,\text {conf}}S^1_{1,in}$, which for the partial-wave amplitudes implies

$$\begin{aligned} t^1_1(s)=t^1_{1,\text {conf}}(s)+t^1_{1,in}(s) +2i\sigma (s)t^1_{1,\text {conf}}\,(s)t^1_{1,in}(s). \end{aligned}$$

(20)

The elastic region is dominated by the $\rho (770)$ resonance and its peak mass will be imposed with an explicit factor in a purely elastic contribution, $t^1_{1,\text {conf}}(s)$, parameterized with a conformal expansion as follows:

$$\begin{aligned} t^1_{1,\text {conf}}(s)&= \frac{1}{\sigma (s)}\frac{1}{\Phi ^1_1(s)-i}, \quad s<1.4 {\mathrm {\,GeV}} \end{aligned}$$

(21)

$$\begin{aligned} \Phi ^1_1(s)&=\frac{\sqrt{s}}{2q^3(s)}(m_\rho ^2-s)\left( \frac{2 m_\pi ^3}{m_\rho ^2 \sqrt{s}}+\sum _{n=0}^N{B_n \omega (s)^n}\right) . \end{aligned}$$

(22)

Once again, $\Phi _1^1(s)=\cot \delta ^1_1(s)$ in the elastic region $s\le 4m_K^2$. As with the S0 wave, the term $\sim 1/\sqrt{s}$ within the parenthesis removes spurious ghosts but makes an almost irrelevant contribution to the fit. For this wave it will be enough to set $N=4$ to obtain a good overall $\chi ^2/d.o.f$ in this region. As before, the conformal variable is defined as

$$\begin{aligned} \omega (s)=\frac{\sqrt{s}-\alpha \sqrt{s_0-s}}{\sqrt{s}+\alpha \sqrt{s_0-s}}, \end{aligned}$$

(23)

but now $s_0=\left( 1.43 {\mathrm {\,GeV}} \right) ^2$ and in order to get an error band whose shape is closer to the actual spread of data, $\alpha$ is chosen so that the expansion center is near the $\pi \pi$ threshold. Values ranging from 0.2 to 0.5 make a suitable parameterization and we use $\alpha =0.3$.

We have already commented that, in contrast to the large $f_0(980)$ effects in the S0 wave, for the P-wave inelastic effects are very tiny below 1.12 ${\mathrm {\,GeV}}$ and very small below $\sqrt{s_m}\equiv 1.4 {\mathrm {\,GeV}}$. Actually, if one is not interested in very high accuracy, using the conformal part of the parameterization alone is almost indistinguishable from using our full partial-wave below 1.12 GeV. However, this very small inelasticity is relevant for the accurate fulfillment of dispersion relations (particularly the $\pi ^+\pi ^0$ Forward Dispersion Relation, see below). Hence, we will include an inelasticity right from $K\bar{K}$ threshold by means of another $t^1_{1, in }(s)$ amplitude defined as

$$\begin{aligned} t^1_{1,in}(s)= & {} \frac{e^{i \delta ^1_{1,in}(s)}-1}{2 i\sigma (s)}, \\ \delta ^1_{1,in}(s)= & {} \bar{J}(s,m_K) \left( K_0+K_1 \frac{s}{m_K^2}\right) \frac{q_\pi ^3(s)}{\sqrt{s} \, m_\pi ^2} \frac{q_K^2(s)}{m_K^2}, \nonumber \end{aligned}$$

(24)

where $q_i=\sqrt{s/4-m_i^2}$ are the CM momenta of the two-pion or two-kaon system and the analytic function $\bar{J}(s,m_i)$ was defined in Eq. (9). Thus, the whole $t^1_1$ amplitude in Eq. (20) has the appropriate kinematic behavior. Namely, the elasticity behaves as $q_K^3(s)/\sqrt{s}$ near $K \bar{K}$ threshold, whereas the phase shift behaves as $q_\pi ^3(s)/\sqrt{s}$ near the $\pi \pi$ one. In addition, since $\bar{J}(s,m_K)$ is real below $K\bar{K}$ threshold, this ensures that the whole $t^1_1$ is elastic for $s<4 m_K^2$ (no inelasticicity has been observed there). Above $K\bar{K}$ threshold $\bar{J}(s,m_K)$ has an imaginary part and therefore both $t^1_{1, in}$ and the whole $t^1_1$ become inelastic. The advantage of using the $\bar{J}$ function is that this parameterization is analytic in the whole energy region from $s=0$ up to 1.4 GeV and provides a straightforward analytic continuation to the complex plane that the usual step functions do not provide. Thus we can also continue $t^1_1$ to the complex plane within its Lehmann ellipse. Note also that $t^1_{1, in}$ contributes with a tiny phase shift below the $K \bar{K}$ threshold, which is on average more than two orders of magnitude smaller than the one coming from the conformal mapping. Thus, as we commented above, in the elastic region $t^1_{1, \text {conf}}$ by itself alone gives a remarkably good description of the whole partial wave. However, in order to get to 1.4 GeV the full Eq. (20) is needed.

In practice we have found that just two constants $K_0, K_1$ together with the conformal parameterization in Eq. (21) are enough to describe the phase shift and inelasticity in the real axis below 1.4 ${\mathrm {\,GeV}}$ as well as the complex plane of the P-wave within the Lehmann ellipse, including the $\rho (770)$ pole obtained in the [28] dispersive analysis.

3.2.2 P-wave parameterization above 1.4 ${\mathrm {\,GeV}}$

As before with the S0 wave, above $\sqrt{s_m}\equiv 1.4 {\mathrm {\,GeV}}$ we will provide just phenomenological fits to the P-wave data, ensuring a continuous matching for the phase and elasticity as well as their derivatives. The matching procedure is similar to that for the S0 wave.

For the P-wave the phase shift will be described using Chebyshev polynomials again. Once the matching with the previous parameterization below 1.4 GeV is implemented, we can obtain $\delta ^1_1(s_m)$ and $\delta ^{\prime \,1}_1(s_m)$, so that the phase shift in the region above 1.4 GeV reads:

$$\begin{aligned} \delta ^1_1(s)&=\delta ^1_1(s_m)+\Delta ^1_1 [x_1(\omega _2(s))+1]+d_0 [x_2(\omega _2(s))-1] \nonumber \\&\quad d_1 [x_3(\omega _2(s))+1], \nonumber \\ \Delta ^1_1&=\frac{\delta ^{\prime \,1}_1(s_m)}{\omega '_2(s_m)}-d_0 x'_2(-1) -d_1 x'_3(-1) \nonumber \\&= \frac{\delta ^{\prime \,1}_1(s_m)}{\omega '_2(s_m)}+4d_0-9d_1. \end{aligned}$$

(25)

We have considered up to three Chebyshev polynomials because, when fitting the data on the P-wave phase above 1.4 ${\mathrm {\,GeV}}$ in the next section, we will just need two degrees of freedom, $d_0$ and $d_1$, to obtain an acceptable $\chi ^2/d.o.f$ for all solutions. The polynomials variable is the same as in the S0 case:

$$\begin{aligned} \omega _2(s)&=2\frac{\sqrt{s}-\sqrt{s_m}}{2 {\mathrm {\,GeV}}-\sqrt{s_m}}-1. \end{aligned}$$

(26)

Also, as it happened in the S0 case, for the calculation of uncertainties we will keep $\delta '^1_1(s_m)$ fixed to its central value.

As we did for the S-wave above 1.4 GeV, for the P-wave elasticity we will use an exponential with a negative exponent to ensure $0\le \eta ^1_1\le 1$. This time we have found that using up to the fourth Chebyshev polynomial is good enough to describe this exponent in all cases and do not produce unwanted oscillations. Thus we write:

$$\begin{aligned} \eta ^1_1(s)=\exp \left[ -\left( \epsilon _0+\sum _{k=1}^4 \epsilon _k \Big (x_k(\omega _2(s))-(-1)^k\Big )\right) ^2\right] . \end{aligned}$$

(27)

Once again, continuity at the matching point fixes

$$\begin{aligned} \epsilon _0=\sqrt{-\log (\eta ^1_1(s_m))}, \end{aligned}$$

(28)

whereas the continuity of the derivative imposes

$$\begin{aligned} \epsilon _1&=-\frac{\eta ^{\prime \,1}_1(s_m)}{2 \epsilon _0 \eta ^1_1(s_m)\omega '_2(s_m)}-\epsilon _2 x'_2(-1)-\epsilon _3 x'_3(-1)\nonumber \\&\quad \,\, -\epsilon _4 x'_4(-1)\nonumber \\&=-\frac{\eta ^{\prime \,1}_1(s_m)}{2 \epsilon _0 \eta ^1_1(s_m)\omega '_2(s_m)}+4\epsilon _2-9\epsilon _3+16\epsilon _4. \end{aligned}$$

(29)

Thus, in practice, three free parameters are needed at most. Once again we remark that the logarithm in Eq. (28) appears in the constants needed for the smooth matching but it does not introduce any spurious analytic structure. For the calculation of uncertainties we will keep ${\eta ^{\prime 1}_1}(s_m)$ fixed to its central value. Thus the central value of the derivative is continuous but its uncertainties might show a small kink.

4 Determination of parameters

The aim of this work is to provide a relatively simple global description for each one of the S0 and P waves of $\pi \pi \rightarrow \pi \pi$ scattering, incorporating all analytic constrains at low energies, including Adler zeros, while also describing the existing data up to 2 GeV. They should also be consistent with the dispersive analysis of data up to 1.4 GeV in [28]. Moreover, such parameterizations should provide also simple but realistic estimates of the uncertainties. In the previous section we have provided such simple and ready to use parameterizations. In this sections we will determine the value of their parameters.

Let us recall that the CFD parameterizations of $\pi \pi \rightarrow \pi \pi$ scattering partial waves obtained in [28] are data fits constrained to fulfill a group of forward dispersion relations up to 1.43 ${\mathrm {\,GeV}}$, together with the more sophisticated Roy and GKPY equations for the partial waves, applicable up to roughly 1.1 ${\mathrm {\,GeV}}$. However they were parameterized with piecewise functions and we now want to mimic them and their uncertainties with a global parameterization. Thus, in the real axis for $4m_\pi ^2<s<s_m$ the CFD partial waves of [28] will be fitted. The CFD has much smaller uncertainties than the output of the dispersion relations themselves, and this is why it is preferred to build a more accurate result. We will impose just the phase shift up to the inelastic $K \bar{K}$ threshold, and both the phase shift and elasticity above it.

In addition, we want our new parameterization to be consistent with the dispersive result in the subthreshold region and in the complex plane, particularly with the resonance pole positions and residues. The CFD are piecewise functions and although some of the pieces contain fair approximations to the poles, they do not provide accurate results in the complex plane. Therefore below the elastic threshold, and in the complex plane, we will fit our global parameterization to the output of GKPY equations, which produces narrower errors than that of Roy equations, while both are compatible among themselves in the whole complex plane and real axis. The fit will run from about $\mathrm {Re}\,s \sim (0 {\mathrm {\,GeV}})^2$ to $\mathrm {Re}\,s \sim (1.12 {\mathrm {\,GeV}})^2$, but always inside the applicability region of the GKPY or Roy dispersion relations, which can be found in [24, 42, 49]. Using such a vast region we are able to describe the scattering lengths, the Adler zeros in the S0 wave, and the $\sigma /f_0(500)$ and $\rho (770)$ pole positions and couplings (the $f_0(980)$ is fixed as input). Due to the smaller uncertainties in the real axis, the final errors bars of our parameterization in the complex plane are smaller than the dispersive ones.

All these features will be imposed on our parameterization by means of a $\chi ^2/d.o.f.$ function, over a grid of points separated by 10 MeV both in the real and imaginary directions within the GKPY/Roy equations applicability region. The input values and uncertainties in this $\chi ^2$ are those of the CFD in the physical region below $s_m$ and of the GKPY output [28] in the subthreshold region and in the rest of the complex plane. Nevertheless, the statistical meaning of the $\chi ^2/d.o.f.\sim 1$ loses part of its purpose, as the results coming from dispersion relations are smooth functions instead of normally distributed points, and their uncertainties are totally correlated between bins. As a result, a value lower than 1 is frequently expected, and we will consider all results below or around 1 as good descriptions of our dispersion relations.

Finally, let us recall that above 1.43 ${\mathrm {\,GeV}}$ no dispersive result exists, thus we will make use of the available experimental data. As explained in Sect. 2 the data sources in this energy region produce three different plausible solutions: the first one, called solution I in this work will fit data from [2, 3, 5, 64]. The second one, that we will call solution II, fits data from a later reanalysis by the CERN-Munich collaboration [4]. Finally, solution III uses a recent update [58] of the (− + −) data solution in [4].

4.1 S0-wave fit

We show in Fig. 1 solutions I, II and III of our new S0-wave parameterization up to 1.9 GeV. Their parameters are listed in Table 1 for solution I, Table 2 for solution II and Table 3 for solution III. By construction, they are almost identical up to 1.4 GeV. Nevertheless, there is an almost imperceptible deviation between them in the inelastic region below 1.4 ${\mathrm {\,GeV}}$ due to their matching to different solutions above 1.4 GeV. Actually, above this splitting point the solutions are fairly different either on their phase, elasticity or both.

It is worth noticing that the uncertainties of solution II are larger for the phase shift, due to the scarcity of data above 1.5 ${\mathrm {\,GeV}}$. In addition, the elasticity data forces solution II to drop first and then to raise in the region between 1.5 and 1.8 GeV, which is hard to explain in terms of known resonances. Above 1.8 GeV there are no data for this solution and our functional forms would give even more oscillations, for which there is no evidence. Thus, to avoid further oscillatory behavior above 1.8 GeV, we have included four elasticity data points above that energy coming from [2, 3] which have huge uncertainties but stabilize the fit. In contrast solution I slowly becomes more and more inelastic as the energy increases which is more natural if more and more channels are open. The phase motion of solution III and the relatively sharp dip of the elasticity, which are quite different from the other solutions, are hints of the presence of the $f_0(1500)$ resonance.

Fig. 1

Comparison of solutions I, II and III (Tables 1, 2, 3) versus data. The gray, blue and green bands correspond to the uncertainty of solutions I, II and III, respectively. Above 1.4 GeV, solution I fits the data of [5, 64] (solid circles) and [2, 3] (solid squares), solution II fits [4] (solid diamonds) and solution III fits the updated (− + −) data from [58] (hollow diamonds). The data coming from [9] (empty squares) and [65] (empty circles) for the phase shift and [66] (solid triangle up),  [67](solid triangle down), [6] (empty squares), [65] (empty circles), [68] (empty triangle up) and [69] (empty triangle down) for the elasticity are just shown for comparison. The red-dashed vertical line separates the region where the fits describe both data and dispersion relation results, from the region above, where the parameterization is just fitted to data. The blue-dotted vertical line stands at the energy of the last data point of solutions II and III

Table 1 Fit parameters of the global parameterization for the S0-wave solution I. $s_p$ is the $f_0(980)$ pole position from the dispersive analysis in [36]

Table 2 Fit parameters of the global parameterization for the S0-wave solution II. $s_p$ is the $f_0(980)$ pole position from the dispersive analysis in [36]

Table 3 Fit parameters of the global parameterization for the S0-wave solution III. $s_p$ is the $f_0(980)$ pole position from the dispersive analysis in [36]

Fig. 2

Comparison between the CFD fit in [28] (blue) and solution I (Table 1, orange band). The energy region dominated by the $f_0(980)$ pole is delimited between the red dashed lines

Concerning the compatibility with the dispersive results in [28], we show in Fig. 2 the comparison between the CFD analysis of [28] and our solution I. Up to 1.4 GeV it is enough to refer to solution I as the global solution, because it is the simplest and all them are almost indistinguishable below 1.4 GeV. The relevant observation from Fig. 2 is that the piecewise CFD and our new parameterization look almost the same below the $K \bar{K}$ threshold and are also very similar and compatible above it. The sharp structure in the region between the two vertical lines in Fig. 2 is dominated by the $f_0(980)$ contribution that we have factored out explicitly in our global parameterization.

All in all, this new parameterization is consistent with the GKPY dispersive data analysis, its output in the complex plane, as well as with the threshold parameters, the Adler zero, the positions of both $\sigma /f_0(500)$ and $f_0(980)$ poles, and the inelastic region up to 1.43 ${\mathrm {\,GeV}}$, which was consistent with Forward Dispersion Relations. This consistency is illustrated in Table 4 where we show the $\chi ^2/d.o.f.\equiv \hat{\chi }^2$ of our fit with the new parameterization in different regions: $\hat{\chi }_1^2$ from $\pi \pi$ to $K\bar{K}$ threshold, $\hat{\chi }_2^2$ from $K\bar{K}$ threshold to 1.4 GeV, $\hat{\chi }_{\mathbb {C}}^2$ in the complex plane within the applicability region, $\hat{\chi }^2_{\delta }$ for the phase above 1.4 GeV and $\hat{\chi }^2_{\eta }$ for the elasticity above 1.4 GeV. All of them are smaller or equal to one for any of our three solutions.

Moreover, this S0-wave global parameterization up to 1.4 GeV is fully consistent with the dispersion relations described in the the GKPY dispersive analysis  [28]. Of course, for such calculation we also need the input for other partial waves and high-energy input described in  [28] and thus we have relegated this discussion to appendix A.

4.1.1 Poles, couplings and low energy parameters

As explained above, the global parameterization is also constrained to describe the dispersive results in the whole complex energy-squared plane. This produces a stable and accurate description of the $\sigma /f_0(500)$ resonance parameters. Actually, in Fig. 3 we show our parameterization and its uncertainties in the first Riemann sheet of the complex plane, which reproduces the output of GKPY equations. In order to see the consistency with the GKPY dispersive result, in the upper panel of Fig. 4 we show the absolute values of the differences between the real part of our new parameterization and the GKPY result divided by the uncertainty of the latter. In the lower panel we show a similar plot for the imaginary parts. Note that our new parameterization lies within the uncertainties of the GKPY for the most part of the region. The only place where there are sizable differences beyond two standard deviations is for Im$\,t_0^0$ in the real axis around 0.9 GeV, but this is the matching point of the two pieces of the CFD parameterization, whereas the GKPY output is much smoother. Thus, our two inputs are slightly incompatible around that region and our new parameterization lies somewhere between both of them.

Fig. 3

Real (top) and imaginary (bottom) parts of the scalar-isoscalar partial wave in the first Riemann sheet of the complex-s plane, within the applicability region of GKPY/Roy equations. There are actually three surfaces on each plot: one for the central value, one for the upper uncertainty and another one for the lower uncertainty band. Note that the behavior of the parameterization is smooth and the uncertainties are small compared to the typical scale of the analytic structures, even deep in the complex plane. We plot solution I, since solutions I, II and III are almost identical in this region

Fig. 4

Within the applicability region of GKPY equations in the complex-s plane, we show the absolute value of the differences between the real (top) and imaginary (bottom) parts of the global parameterization and the GKPY equations, divided by the uncertainty of the latter. We plot results for solution I but the other two are identical in this region

In addition, we list in Table 5 the parameters of both the $\sigma /f_0(500)$ and $f_0(980)$ resonances compared to their GKPY dispersive values in [28]. It is worth noticing that the uncertainties of the $\sigma /f_0(500)$ resonance associated to this fit are a bit smaller than the GKPY determination [36]. This is due to the fact that besides the GKPY output in the complex plane, we are fitting the CFD in the real axis, which has smaller uncertainties.

As for the $f_0(980)$ resonance, we have included in the fit the pole position obtained by means of the GKPY equations in [36]. The main reason is that phenomenological fits cannot extract its accurate parameters in a very stable way (see for instance [70]). In particular, the CFD fit of [28] does not provide an accurate estimate of its position and one has to rely on the numerical dispersive approach. However, with our new parameterization, the $f_0(980)$ is no longer a problem, as both the data, the cusp effect and the pole position are factored out into a simple, yet versatile functional form. Once again, the coupling of the $f_0(980)$ to $\pi \pi$ has smaller uncertainties than the GKPY determination (Table 5), since the CFD partial wave, with its small uncertainty, is also fitted in the real axis to obtain our new parameterization.

Table 4 Results in terms of $\hat{\chi }^2$ of the S0 solutions I, II and III in different regions

Last, but not the least, the global parameterization yields relatively accurate threshold and sub-threshold parameters (like the Adler zero), which are included in Table 6. These values are compatible with those of the dispersive data analysis of the Madrid-Krakow group [28] and therefore also with the dispersive analysis matched to two-loop ChPT of the Bern group [24, 25].

4.2 P-wave fit

Following the same procedure just applied to the S0 wave, in the physical region and below 1.4 ${\mathrm {\,GeV}}$ we will fit our P-wave parameterization to the CFD P-wave of [28]. Note that this P-wave parameterization describes data from both $\pi \pi$ scattering [4, 6, 71] and the pion vector form factor from [72, 73], while fulfilling at the same time the GKPY/Roy equations up to 1.12 GeV and Forward Dispersion Relations up to 1.43 GeV.

Table 5 Pole positions and $\pi \pi$ couplings of both $f_0(500)$ and $f_0(980)$ resonances from our global parameterization. Almost indistinguishable values would be obtained for solutions I, II and III. Note that they are very compatible with the GKPY dispersive results in [36]

Table 6 Adler zero and threshold parameters. The latter in customary $m_\pi =1$ units. They are almost indistinguishable for solutions I, II and III

While in recent years there have been no new $\pi \pi$-scattering experimental analyses, the pion-vector form factor data described in [28] is nowadays outdated with respect to the modern and precise measurements in [74, 75]. These experimental data have been recently considered as input in a dispersive analysis for $e^+e^-\rightarrow \pi ^+\pi ^-$ [76], where the $\pi \pi$-scattering P-wave phase shift is obtained from the Roy-equation analysis in [25, 31] and the P-wave phase-shift values at $\sqrt{s}=0.8$ and $1.15 {\mathrm {\,GeV}}$ are considered as fit parameters. In this way, the result depicted in Fig. 16 in [76] should be taken as the most precise and updated dispersive determination of the $\pi \pi$ scattering P-wave phase-shift. Nevertheless, this recent determination is compatible with the CFD result in [28] within uncertainties¹ and, for consistency, it will be still used as our input in the real axis.

Once more, in the subthreshold region and in the complex plane we will fit the GKPY-equation dispersive results. As done for the scalar channel, we will only consider the energy region within the Lehmann ellipse, where both Roy and GKPY equations are formally valid. Above 1.43 ${\mathrm {\,GeV}}$ there are no further dispersive results and hence we will only describe the available experimental data, which come from a single scattering experiment performed by the CERN-Munich Collaboration. In addition, in the vector case there is a relevant difference between the best solution of the original CERN-Munich result published in 1973 [2] and the (− − −) solution of the 1975 collaboration reanalysis [4]. The revisited and modified (− + −) solution in [58] lies somewhat in between.

The behavior of the original P-wave result shows a large interference in the region between 1.5 and 1.8 ${\mathrm {\,GeV}}$. Namely, within these 300 MeV, the phase shift changes by more than $20^\circ$ and the elasticity, starting from almost 1, decreases to less than 0.5 to return back to 1. This behavior could only be explained if the $\rho '$ and $\rho ''$ resonances and the $K \bar{K}$ channel would interfere strongly, which is in contradiction with the experimental values for the width and couplings of these two resonances [77, 78]. Thus, the solution (− − −) of Hyams 75 [4] is the one customarily used in the literature. However, we will fit three solutions for completeness, as we have done for the S0-wave. The original CERN-Munich result [2, 3] will be called solution I, whereas the fit to the updated reanalysis of [4] will be called solution II and the fit to the updated (− + −) solution in [58] will be called solution III.

Table 7 Results in terms of $\hat{\chi }^2$ of the P-wave solutions I, II and III, in different regions

As previously done for the S0 wave we will fit our P-wave global parameterization described in Sect. 3.2 to a 10 MeV-spaced grid of GKPY output values within their applicability region in the complex plane and to the CFD parameterization in the real axis at energy points separated by 5 MeV. In addition we add the $\chi ^2/d.o.f.$ of the data above 1.4 GeV, although for the phase shift of solutions I and II we have added $1^\mathrm{o}$ as a systematic uncertainty, since the nominal uncertainties in some regions are unrealistically small, particularly for solution II. The fit minimizes a $\chi ^2/d.o.f.$ function whose uncertainties are those of the GKPY or the CFD partial wave. Once more, even though our $\chi ^2/d.o.f.$ does not have a well-defined statistical meaning, it ensures a nice description of the input as seen in the $\hat{\chi }^2\equiv \chi ^2/d.o.f.$ values, given in Table 7. They come out $\hat{\chi }^2\simeq 1$ or less in all regions (we follow the same notation as for the S0 wave).

Fig. 5

Comparison between our three P-wave solutions and scattering data. The gray, blue and green bands correspond to solution I, II and III, respectively. The red dashed vertical line separates the region where the fits describe both data and dispersive results, from the region above where the parameterization is just fitted to the data. Namely, above 1.4 GeV solution I is fitted to  [2, 3] (solid diamonds), solution II to [4] (solid upward triangles) and solution III to [58] (solid downward triangles). The data from  [6] (solid squares),  [71] (solid circles) are just shown for completeness. The blue dotted vertical line depicts the energy of the last data point of solutions II and III

The resulting P-wave phase shift and elasticity are plotted in Fig. 5 and their parameters are collected in Tables 8, 9 and 10 for solutions I, II and III, respectively. Both from the figure and the tables we can see that they are almost identical up to 1.4 GeV. The uncertainties are described by the gray band for solution I, the blue band for solution II, and the green band for solution III. Below $s_m$ the three of them reproduce nicely the uncertainties given in [28].

Table 8 Fit parameters of the Global parameterization for the P-wave solution I

Table 9 Fit parameters of the Global parameterization for the P-wave solution II

Concerning the region below 1.4 GeV, note that the CFD uncertainties are extremely small below $K\bar{K}$ threshold. For this reason, in order to ensure an accurate description of the error band in this region, we chose $\alpha =0.3$ in Eq. (21) so that the center of the conformal expansion in our new amplitude is close to the $\pi \pi \rightarrow \pi \pi$ threshold. In this way, the uncertainties there are dominated by the lowest conformal parameters $B_0$ and $B_1$, ensuring that the values of the threshold parameters, given in Table 11, are also consistent with the dispersive values in [28]. In contrast, the new parameterization uncertainties close to 1.4 ${\mathrm {\,GeV}}$ are smaller than those of the CFD in [28], which is a consequence of describing simultaneously the experimental data up to 2 ${\mathrm {\,GeV}}$. The P-wave elasticity in [28] is compatible with 1 below 1.12 ${\mathrm {\,GeV}}$ and slightly smaller below 1.4 ${\mathrm {\,GeV}}$. This is why it can be reproduced by using only two free constants in the parameterization given in Eq. (24).

The $\rho (770)$-pole parameters are given in Table 12 and are identical for all solutions. Central values and uncertainties are nicely compatible with the dispersive results in [36].

Furthermore, as we already commented for the S-wave, this P-wave global parameterization up to 1.4 GeV is fully consistent with all the dispersion relations described in the the GKPY dispersive analysis  [28]. Since such a calculation requires as input the other partial waves and high-energy input described in  [28], we have relegated it to Appendix A.

Above the matching point at 1.4 ${\mathrm {\,GeV}}$ the three solutions coming from the CERN-Munich experiment are incompatible among themselves. The behavior of solution I suggests a strong interference between the $\rho '$ and $\rho ''$, with a sizable phase change around 1.6 ${\mathrm {\,GeV}}$ and a dip structure in the elasticity at the same energy. In contrast, solutions II and III look smoother. Namely, the phase grows slowly above $180^{^{\circ }}$ and the elasticity has a less pronounced dip. In addition, the uncertainties quoted in [4] and [58] are slightly smaller, which leads to more constrained uncertainty bands. Nevertheless we emphasize once more that above 1.4 GeV, we consider our parameterizations purely phenomenological.

Table 10 Fit parameters of the global parameterization for the P-wave solution III

Table 11 P-wave threshold parameters in customary $m_\pi =1$ units. They are almost indistinguishable for solutions I, II and III

Table 12 Pole position and $\pi \pi$ coupling of the $\rho (770)$, which are almost indistinguishable for solutions I, II and III. They nicely agree with the GKPY dispersive result in [36] that we also provide for comparison

5 Summary

In this work we have provided a global parameterization of the data for each one of the S0 and P-waves of $\pi \pi \rightarrow \pi \pi$ scattering up to almost 2 GeV. We have made an explicit effort to keep it relatively simple in order to be easy to implement in further phenomenological and experimental analyses (in final state interactions, isobar models, etc...).

The advantages of these parameterizations are that they describe experimental data up to 2 GeV consistently with the dispersive representation in [28] and its uncertainties up to its maximum applicability region of 1.4 GeV. In addition, they reproduce the dispersive results within their applicability region in the complex-s plane, as obtained in [36], including the poles associated to the $\sigma /f_0(500)$, $f_0(980)$ and $\rho (770)$ resonances. Moreover, their low-energy behavior is compatible with the dispersive results for the threshold parameters and the S0 Adler zero and hence with the constraints due to the QCD spontaneous chiral symmetry breaking.

Actually, these new parameterizations reproduce the results and uncertainties of a previous piecewise fit that was constrained to satisfy Forward Dispersion relations up to 1.4 GeV and partial-wave dispersion relations (Roy and GKPY equations) up to 1.12 GeV. The latter were used in [28] to obtain a rigorous analytic continuation to the complex plane which, together with its uncertainties, is also described when continuing analytically our new parameterization, without the need for a numerical integration of the dispersion relations. This is why the pole positions and residue of the $\sigma /f_0(500)$, the $f_0(980)$ and the $\rho (770)$ are so well implemented. It also allows our parameterization to be used consistently in applications with isobar models, so popular in experimental analyses.

The new parameterizations also reproduce the existing data from 1.43 to 2 GeV, although the dispersion relations do not reach these energies. Moreover, in this region, there are three contradictory data sets, and we thus provide three solutions for each wave that describe phenomenologically either one of the conflicting sets. Nevertheless, below 1.4 GeV these three solutions agree to the point of being almost indistinguishable and are consistent with the dispersive analysis.

We hope that the simplicity and the remarkable analytic properties of this data parameterization can be useful for future phenomenological and experimental studies in which $\pi \pi \rightarrow \pi \pi$ interactions are needed.

Appendix A: Roy-Steiner and forward dispersion relations

As we have explained in the main text, the global parameterizations we present in this work mimic very well the dispersive piece-wise CFD parameterization of [28] and their uncertainties up to almost 1.4 GeV. Of course, they are not exactly identical to the CFD and one might wonder if the dispersion relations used in [28] to constrain the piece-wise parameterization are still satisfied. In this appendix we show that this is indeed the case.

Here we only show the results for solution I since, as we have emphasized repeatedly, up to 1.4 GeV all our three solutions are almost indistinguishable. The explicit equations for each dispersion relation as well as the rest of the input apart from the S0 or P-waves below 1.4 GeV can be found in [28] and references therein. For the S0 and P waves below 1.4 GeV we use the expressions in the main text here. Since our DR only reach up to 1.4 GeV at most, the results with the other solutions are indistinguishable.

Let us simply recall that each dispersion relation provides the real part of the amplitude as an integral from threshold to infinity of the imaginary part times some specific kernel for each relation. Sometimes some subtraction terms (polynomials in s), which are real, are added to this integral, and their coefficients (subtraction constants) are obtained from the values of the amplitude at fixed values, typically at or around threshold. The imaginary part is taken from the fits to data. That integral and the subtraction terms provide what we call the “Dispersive” result. When it agrees within uncertainties with the values of the amplitude obtained from our fits, we consider that the dispersion relation is well satisfied.

Thus, in Fig. 6 we show the “Dispersive” result versus our fit for the the S0, P and S2 Roy equations. It can be seen that in all cases they are well satisfied within uncertainties within the whole applicability region up to $\sqrt{s}\sim 1.1\,$GeV. The same can be said about the GKPY equations, which are similar to Roy equations but with just one subtraction. They are more sensitive to the resonance region and less to the threshold region. In any case, it can be observed in Fig. 7 that they are also well satisfied within uncertainties within their applicability region.

Finally, we show in Fig. 8 the fulfillment of Forward Dispersion Relations (FDR). It should be noted that these are not written for partial waves, but for the amplitudes

$$\begin{aligned} F^{0}_{0}= & {} \frac{1}{3}(F^{(0)}+2F^{(2)}),\quad F^{0}_{+}=\frac{1}{2}(F^{(1)}+F^{(2)}), \\ F^{I_t=1}= & {} \frac{1}{3}F^{(0)}+\frac{1}{2}F^{(1)}-\frac{5}{6}F^{(2)}. \end{aligned}$$

where

$$\begin{aligned} F^{(I)}(s,t)=\frac{8}{ \pi } \sum _\ell (2\ell +1)P_\ell (\cos \theta ) t^{I}_\ell (s). \end{aligned}$$

(30)

Once again we see that these FDR are well satisfied up to 1.4 GeV, which is the energy up to which they were studied before in [28]. Above that energy, Regge theory is being used as input for the integrals instead of partial waves. Note that in [28] the FDRs were checked up to 1.42 GeV whereas here we only check them up to 1.4 GeV since it is at that point that we make the matching with the different solutions of higher energies and our new global parameterization starts deviating from our piece-wise solution to accommodate the different data sets above 1.4 GeV. As seen in Fig. 8 the agreement is remarkable except in the final points very close to 1.4 GeV, i.e, the matching point, where our global fit has to deviate slightly from the old CFD to match the data above 1.4 GeV. Nevertheless, note that, even in [28] there was a small deviation in the FDR around 1.4 GeV due to the discontinuous matching with the Regge description.

Fig. 6

Results for Roy equations. Blue lines: real part coming from our new parameterizations. Orange lines: the result of the dispersive integrals. The gray bands cover the uncertainties in the difference between both. From top to bottom: a S0 wave, b P wave, c S2 wave

Fig. 7

Results for GKPY equations. Blue lines: real part coming from our new parameterizations. Orange lines: the result of the dispersive integrals. The gray bands cover the uncertainties in the difference between both. From top to bottom: a S0 wave, b P wave, c S2 wave

Fig. 8

Results for forward dispersion relations. Blue lines: real part coming from our new parameterizations. Orange lines: the result of the dispersive integrals. The gray bands cover the uncertainties in the difference between both. From top to bottom: a the $\pi ^0\pi ^0$ FDR, b the $\pi ^0\pi ^+$ FDR, c the FDR for $I_t=1$ scattering

Moreover, following [28] and to quantify the agreement between the fit and the dispersive result we have defined $\Delta _i(s)$ as the difference between the dispersive output and the fit for each dispersion relation i, whose uncertainties are called $\delta \Delta _i(s)$. Next, the average discrepancies are defined as follows:

$$\begin{aligned} \bar{d}_i^2\equiv \frac{1}{\hbox {number of points}} \sum _n\left( \frac{\Delta _i(s_n)}{\delta \Delta _i(s_n)}\right) ^2 , \end{aligned}$$

(31)

where the values of $\sqrt{s_n}$ are taken at intervals of 25 MeV as in [28]. These discrepancies are similar in form to a $\chi ^2$ function and when they are $\sim 1$ or smaller it means that the dispersion relation is well satisfied. The resulting discrepancies are shown in Table 13 for different energy intervals. The fulfillment of these dispersion relations is remarkable, since all the average discrepancies are less than one except for the 1.25 value of the $\pi ^+\pi ^0$ up to 1.4 GeV. Even that is very good taking into account that the uncertainties of the global fit are somewhat smaller close to 1.4 GeV than for the old CFD piece-wise parameterization, due to the matching with additional data above that energy. Actually these discrepancies are on the average very similar in size to those obtained in [28] with the piece-wise CFD parameterizations. The consistency of our global parameterization with respect to dispersion relations is thus remarkable and comparable to the one in [28].

Table 13 Average discrepancies $\bar{d}_i^2$ of our global parameterizations for each dispersion relation

As a final comment, note that if we had neglected the P-wave tiny inelasticity between 1 GeV and 1.12 GeV, setting $\eta ^1_1=1$ there, then $\bar{d}_{\pi ^+\pi ^0}^2$ up to 1.4 GeV would have been $\sim 6$. It is for this reason that we have have chosen to introduce $t^1_{1,in}$ in our P-wave parameterization, Eq. (20).