A new parametrization for the scalar pion form factors

We derive a new parametrization for the scalar pion form factors that allows us to analyze data over a large energy range via the inclusion of resonances, and at the same time to ensure consistency with the high-accuracy dispersive representations available at low energies. As an application the formalism is used to extract resonance properties of excited scalar mesons from data for B¯s0→J/ψππ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\bar{B}}^0_s\rightarrow J/\psi \pi \pi $$\end{document}. In particular we find for the pole positions of f0(1500)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0(1500)$$\end{document} and f0(2020)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_0(2020)$$\end{document}1465±18-i(50±9)MeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1465\pm 18 - i (50\pm 9)\,\text {MeV}$$\end{document} and 1910±50-i(199±40)MeV\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1910\pm 50 - i(199\pm 40)\,\text {MeV}$$\end{document}, respectively. In addition, from their residues we calculate the respective branching ratios into ππ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi \pi $$\end{document} to be (58±31)%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(58\pm 31)\%$$\end{document} and (1.3±1.8)%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(1.3\pm 1.8)\%$$\end{document}.


Introduction
The scalar isoscalar sector of the QCD spectrum up to 2 GeV has been of high theoretical and experimental interest for many years. One of the main motivations for these investigations is the hunt for glueballs: their lightest representatives are predicted to occur in the mass range between 1600 and 1700 MeV with quantum numbers 0 ++ [1][2][3][4]. The most straightforward way to identify glueball candidates is to count states with and without flavor quantum number and see if there are supernumerary isoscalar states; see, e.g., the minireview on non-qq states provided by the Particle Data Group (PDG) [5] or the reviews Refs. [6,7]. Unfortunately, regardless of the year-long efforts, the scalar isoscalar spectrum is still not fully resolved: e.g. there is still an ongoing debate whether the f 0 (1370) exists or not [6]. One problem might be that most analyses of experimental data performed so far are based on fitting sums of Breit-Wigner functions, which can lead to reaction-dependent results. To make further progress, it therefore appears compulsory to employ a e-mail: ropertz@hiskp.uni-bonn.de b e-mail: c.hanhart@fz-juelich.de c e-mail: kubis@hiskp.uni-bonn.de parametrizations that allow one to extract pole parameters, for those by definition do not depend on the production mechanism. This requires amplitudes that are consistent with the general principles of analyticity and unitarity. In this paper we present a new parametrization for the scalar pion form factors that has these features built in, and in addition maps smoothly onto well constrained low-energy amplitudes.
The two-pion system at low energies is well understood from sophisticated investigations based on dispersion theory-in particular the ππ-KK phase shifts and inelasticities can be assumed as known from threshold up to an energy of about s = (1.1 GeV) 2 [8][9][10][11][12][13]. From this information, quantities like the scalar non-strange and strange form factors for both pions and kaons can be constructed, again employing dispersion theory [14][15][16][17][18][19][20][21]. The resulting amplitudes, which capture the physics of the f 0 (500) (or σ ) and the f 0 (980), were already applied successfully to analyze various meson decays, see, e.g., Ref. [20]. In particular the non-Breit-Wigner shape of these low-lying resonances [22] is taken care of automatically. However, to also include higher energies in the analysis, where additional inelastic channels become non-negligible and higher resonances need to be included, one is forced to leave the safe grounds of fully model-independent dispersion theory and to employ a model. Ideally this is done in a way that the amplitudes match smoothly onto those constructed rigorously from dispersion relations. Moreover, to allow for an extraction of resonance properties, the extension needs to be performed in a way consistent with analyticity.
A formalism that has all of these features was introduced for the pion vector form factor in Ref. [23]. In that case, the low-energy ππ interaction can safely be treated as a single-channel problem in the full energy range where high-accuracy phase shifts are available, since the two-kaon contribution to the isovector P-wave inelasticity is very small [12,24]. 1 However, this is not true for the isoscalar S-wave, clearly testified by the presence of the f 0 (980) basically at the KK threshold with a large coupling to this channel [26,27]. Thus, in order to apply the formalism of Ref. [23] to the scalar isoscalar channel it needs to be generalized. This is the main objective of the present article. As an application we test the amplitudes on data forB 0 s → J/ψπ π/KK recently measured with high accuracy at LHCb [28,29], which allows us to extract the strange scalar form factor of pions and kaons up to about 2 GeV and to constrain pole parameters and branching fractions of two of the heavier f 0 resonances in that energy range. This paper is organized as follows. In Sect. 2, we derive the unitary and analytic scalar form factor parametrization to be used. In Sect. 3 we illustrate its application in a coupled-channel analysis of the decaysB 0 s → J/ψπ π and B 0 s → J/ψ KK . Specifically, we discuss the stability of our fits under changing assumptions for the parametrization concerning the number of resonances, the degree of certain polynomials, as well as the approximation in the description of the effective four-pion channel. In addition, in Sect. 4 we extract pole parameters, in particular for both the f 0 (1500) and the f 0 (2020), via the method of Padé approximants for the analytic continuation to the unphysical sheets. The paper ends with a summary and an outlook in Sect. 5.

Formalism
The derivation of the form factor parametrization is presented for the strange scalar isoscalar pion (kaon) form factor Γ s π (Γ s K ). These are related to the matrix elements where s = ( p 1 + p 2 ) 2 . The Γ s i (s), i = π, K , defined this way are invariant under the QCD renormalization group. Since the scalar isoscalar ππ system is strongly coupled to the KK channel via the f 0 (980) resonance, a coupled-channel description becomes inevitable even for energies around s = 1 GeV 2 . In this paper we present a parametrization for the scalar form factors valid at even higher energies. This becomes possible via the explicit inclusion of further inelasticities and resonances. Below 1 GeV the system is strongly constrained by dispersion theory using a coupled-channel treatment of ππ and KK [20]. At higher energies experimental data indicate that further inelasticities are usually accompanied by resonances. We thus derive a parametrization that allows for resonance exchange at higher energies. Those resonances also act as doorways for the coupling of the system to the additional channels. At the same time we make sure that their presence does not distort the amplitude at lower energies. To be concrete, here we consider in addition to ππ (channel 1) and KK (channel 2) an effective 4π channel (channel 3), modeled by either ρρ or σ σ . Three-channel models with an effective σ σ channel have been considered in the literature before [15,30], while some of the f 0 states between 1.3 and 2 GeV have even been hypothesized to be dynamically generated by attractive interactions between ρ mesons [31][32][33][34]. It should become clear from the derivation, however, that the formalism allows for the inclusion of additional channels in a straightforward manner.
The derivation starts from the scalar isoscalar scattering amplitude T (s) i f , where i and f denote the initial-and finalstate channels. To implement unitarity and analyticity we use the Bethe-Salpeter equation, which reads in operator form.
Here V i f denotes the scattering kernel of the initial channel i into the final channel f . The loop operator G is diagonal in channel space and provides the free propagation of the particles of channel m. For example, at the one-loop level the above equation generates an expression of the form for ππ rescattering, with P being the total 4-momentum of the system such that P 2 = s. For m = 1, 2, the discontinuity of the loop operator element G mm reads where σ m (s) = 1 − 4M 2 m /s is the two-body phase space in the given channel, and M m denotes the pion and kaon masses for channels 1 and 2, respectively. For the third channel, we need to include the finite width of the two intermediate (ρ and σ ) mesons; we write where λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + ac + bc) is the Källén function. Here the spectral density for the state k, ρ k (m 2 ), is given as with the energy-dependent width where Γ k (m k ) denotes the nominal width (mass) of the resonance and L k the angular momentum of the decay with L k = 1 and 0 for the ρ and the σ , respectively. The F (L) R (s) denote barrier factors that prevent the width from growing continuously. We employ the parametrization of Refs. [35,36], where their explicit forms for L = 0, 1, 2 are given by , and the hadronic scale r R = 1.5 GeV −1 . Note that as long as no exclusive data are employed for the 4π final state, the amplitudes are not very sensitive to the details how, e.g., the spectral density of the σ meson is parametrized, since it enters only as the integrand in the self energies of the resonances. However, the analysis is somewhat sensitive to the differences between a ρρ and a σ σ self energy, since the energy dependence of the two is quite different, given the different resonance parameters and the different threshold behavior. We come back to this discussion later in this section.
To proceed with the derivation we split the scattering kernel into two parts, V = V 0 + V R , conceptually following the so-called two-potential formalism [37]. The effect of V 0 will eventually be absorbed into the dispersive piece fixed by the low-energy ππ-KK T -matrix input. Its explicit form is needed at no point; one may think of it as the driving term of a Bethe-Salpeter equation Since T 0 it is the solution of a scattering equation, T 0 is unitary. In particular, we may write where δ 0 is the scalar isoscalar ππ phase shift, ψ 0 the phase of the ππ → KK scattering amplitude, and g 0 its absolute value. The inelasticity η 0 is related to g 0 via The effects of resonances heavier than the f 0 (980) enter the amplitude via V R . By means of V R we can construct the resonance T -matrix T R , related to the full T -matrix via T = T R +T 0 . Since T 0 is unitary by itself, T R cannot be independent of T 0 in order to respect the Bethe-Salpeter equation (2). Solving for T R we obtain To proceed, we define the vertex function Ω via Its discontinuity is given by which agrees with the discontinuity of the Omnès matrix derived from the scattering T -matrix T 0 [38,39]. Therefore it can be constructed from dispersion theory: Numerical results for Ω i j (s) based on the T -matrix of Ref. [40] are shown in Fig. 1. One observes in particular the signature of the f 0 (500) or σ -meson, i.e. the broad bump in the imaginary part of Ω 11 (s) below 1 GeV, accompanied by a quick variation of the real part, which clearly cannot be parametrized by a Breit-Wigner form. For an earlier discussion about this fact see Ref. [22]. Using T R = Ω t R Ω t and V 0 GΩ = Ω − 1, which follows from inserting Eqs. (9) into (13), one obtains a Bethe-Salpeter equation for t R , Note that Eq. (16) does not depend on V 0 explicitly. It appears only implicitly, since the loop operator G, describing the free propagation of the two-meson states, needs to be replaced by the dressed loop operator (GΩ), which describes the propagation of the two-meson state in the presence of the interaction T 0 , in order to preserve unitarity. The discontinuity of this self-energy matrix Σ = GΩ is given by The discontinuities of the loop functions for the two-body channels and the 4π channel were given in Eqs. (4) and (5), respectively. Equation (17) allows us to write Σ as a oncesubtracted dispersion integral, The resulting self-energy functions Σ i j (s) are displayed in Fig. 2. The subtraction constants can be absorbed in a redefinition of the yet undefined potential V R . Please observe that the component Σ 33 looks very different for the two different model assumptions employed. For example, the self energy from σ σ intermediate states rises very quickly right from the 4π threshold, while the one for ρρ sets in significantly later. This difference reflects that the σ decays into two pions in an S-wave, while the ρ decays in a P-wave. On the other hand, since the discontinuity of G 33 enters in the expression for Σ 33 only as the integrand, this component of the self energy is not very sensitive to the details of the concrete parametrizations employed for the spectral functions. The full solution for the scattering matrix is thus given by In order to obtain a parametrization for the form factor, we adapt the P-vector formalism [41] to the system at hand. The isoscalar scalar form factor Γ s i is written as where M i is an analytic term describing the transition from the source to the channel i. Inserting the parametrization of Eq. (19) we obtain, after some straightforward algebra, As T 0 captures the physics in the ππ and KK channels at energies below 1 GeV including the f 0 (500), the f 0 (980), and the impact of the corresponding left-hand cuts (left-hand cuts in the other channel(s) are neglected by construction), the potential V R should predominantly describe the resonances above 1 GeV. In order to reduce their impact at low energies, we subtract V R at s = 0 and arrive at The bare resonance masses, m r , as well as the bare resonancechannel coupling constants, g r i , are free parameters that need to be determined by a fit to data. The subtraction constants are effectively absorbed into T 0 that by construction captures all physics close to s = 0.
The most general ansatz for M reads where the parameters c i = Γ s i (0) provide the normalizations of the different form factors. Here the isospin Clebsch-Gordan coefficients were absorbed into the definition of the form factors. This means explicitly while for the third channel we absorb these factors into the coupling constants. The bare resonance masses and the corresponding couplings g r i are the same as before. The parameters α r , which quantify the resonance-source couplings, and the slope parameters γ i are additional free parameters. This completely defines the formalism. Clearly, the number of inelastic channels can be extended in a straightforward way, however, for the concrete application studied in the following section, three channels turn out to be sufficient as long as no exclusive data for additional channels become available.

Parametrization of the decay amplitudes
As an example, we now apply the formalism introduced in the previous section to the decaysB 0 s → J/ψ π + π − (K + K − ), analyzing data taken by the LHCb collaboration [28,29]. The dominant tree-level diagram for the corresponding weak transition on the quark level is displayed in Fig. 3.
It has been argued previously [20,42] that the S-wave projection of the appropriate helicity-0 amplitude forB 0 s → J/ψ M 1 M 2 transitions are proportional to the corresponding strange scalar form factors of the light dimeson sys- The hadronization of thess quark pair into π + π − (S-wave dominated) is given by the scalar form factor Γ s π tem M 1 M 2 ; in particular, there are chiral symmetry relations between the different dimeson channels that fix the relative strengths to be equal to those of the matrix elements in Eq. (1) at leading order in a chiral expansion [42]. We conjecture here that the same will still hold true for the inclusion of the effective third (4π ) channel. In this sense, theB 0 s decays allow to test the pion and kaon strange scalar form factors, up to a common overall normalization.
A previous dispersive analysis [20], which considered the ππ-KK coupled-channel system, worked well in the energy region up to 1.05 GeV. However, due to higher resonances and the onset of additional inelasticities the framework could not be applied beyond this energy. Our new parametrization allows us to overcome this limitation, while it guarantees at the same time a smooth matching onto the amplitudes employed in Ref. [20]. The data are provided in terms of angular moments Y 0 L ( √ s), which are given as angular averages of the differential decay rates where Θ is the scattering angle between the momentum of the dipion system in theB 0 s rest frame and the momentum of one of the pions. We express the decay amplitude in terms of the partial-wave-expanded helicity amplitudes H L λ , where L denotes the angular momentum of the pion or kaon pair, and λ = 0, , ⊥ refers to the helicity of the J/ψ. The angular moments are then given as and for the moments of relevance for this work; see Refs. [20,42] for details. In addition to the pion momentum in the dipion rest frame p π introduced earlier, we also use the J/ψ momentum in theB 0 s rest frame, The scalar helicity amplitude H 0 0 can be related to the scalar isoscalar form factor Γ s i as where the normalization factor N absorbs weak coupling constants and the pertinent Wilson coefficients, as well as meson mass factors and decay constants [20,42]. Here i denotes the relevant channel. For the form factors we use the parametrization introduced in Sect. 2. Since the main focus of our analysis lies on the S-waves, we approximate the Pand D-waves as Breit-Wigner functions [43], for L ≥ 1. The free parameters introduced here are the strength h R λ of the resonance R with helicity λ, its phase φ R λ , and a total rescaling factor w L λ for the helicity amplitude R we use z = r 2 R p 2 π with r R = 1.5 GeV −1 as in Eq. (8) [36]. The position as well as width of the corresponding resonance is then included in the Breit-Wigner function with an energy-dependent width Γ R (s) (7). Since the only interference term in the angular moments considered, Eqs. (26) and (27), is the S-D-wave interference in Y 0 2 , our fits are only sensitive to the relative phase motion of H 0 0 and H 2 0 . To reduce the total number of free parameters for all partial waves except the S-wave, we fix the resonance masses m R as well as their respective widths to the central values found in Refs. [28,29]. Furthermore we fix both h R λ as well as φ R λ with λ = , ⊥ to the central values of the LHCb fits. However, since the phase motion of our S-wave will be different from the one of the LHCb parametrization [20], we allow w R λ to vary. For the helicity amplitude H 2 0 we keep both h R 0 as well as φ R 0 flexible. To avoid unnecessary parameters we set w 2 0 = 1. The number of free parameters is discussed in more detail in Sect. 3.2.

Fits to the decay data
In this section we discuss the fit using the form factor parametrization of Eq. (21) to the data measured forB 0 s → J/ψπ + π − [28] andB 0 s → J/ψ K + K − [29], which are presented as angular moments related to the helicity amplitudes via Eqs. (26) and (27). Note that these angular moments have an arbitrary normalization and need to be rescaled to their physical values. The integrated partial width is given by The correctly normalized angular moments, Y 0 L norm , can be obtained from those published, Y 0 L LHCb , by We determine the partial decay rates and the branching ratios [5] B B 0 s → J/ψ π + π − = (2.09 ± 0.23) × 10 −4 , The dispersive approach using the Omnès matrix already captures the physics of the f 0 (500) and f 0 (980) resonances. In order to extend the description further, we use N R additional resonances. As outlined above, the S-wave contains in total up to (N c + N s + 1)N R + 2N c N s parameters, where N c (N s ) denotes the number of channels (sources) included; in this study N s = 1, N c = 3, and N R is either 2 or 3, depending on the fit. The last term in the sum above comes from the non-resonant couplings of the system to the source. The number of those parameters can be reduced from the observation that the normalizations of the pion and the kaon form factors can be fixed to c 1 = 0 and c 2 = 1 [20]. Since the four-pion channel is expected to couple similarly weakly to anss source as the two-pion one (given OZI suppression at s = 0), we also set c 3 = 0. Thus the only free parameter from the constant terms in the sources M i can be absorbed into the overall normalization N introduced in Eq. (28). Below we present fits without (γ i = 0, resulting in 5N R parameters) as well as with linear terms in the production vertex defined in Eq. (23) (γ i = 0, providing three more free constants).
For the decayB 0 s → J/ψ π π the dipion system is in an isoscalar configuration; due to Bose symmetry the pions can therefore only emerge in even partial waves. Since we restrict ourselves to a precision analysis of the S-wave, we adopt the D-waves of Ref. [28] and accordingly include two resonances, namely f 2 (1270) and f 2 (1525). For the 0 polarization we introduce four new parameters given by the amplitude h R 0 and φ R 0 , while we fix w 0 0 = 1. For the other two helicity amplitudes we constrain h R λ and φ R λ while keeping w 0 λ variable. This gives another two free parameters. In total we obtain six additional free parameters.
Since K + and K − do not belong to the same isospin multiplet, they do not follow the Bose symmetry restrictions. Thus the P-wave in the decayB 0 s → J/ψ K + K − is nonnegligible and, in fact, dominant. It shows large contributions of the φ(1020) as well as of the φ(1680). Since the P-wave does not interfere with S-or D-waves in the angular moments Y 0 0 and Y 0 2 , we adopt the parameters of LHCb [29]. In order to allow for some flexibility, we also fit w 1 λ , resulting in three parameters. The D-wave includes the resonances f 2 (1270), f 2 (1525), f 2 (1750), and f 2 (1950). For λ = 0 we fit both h R 0 as well as φ R 0 with fixed w 2 0 = 1, resulting in eight free parameters. For the other helicity amplitudes we stick to the LHCb parametrization and keep w 2 λ free, which results in two additional fit parameters. Therefore in total we have 13 additional free parameters for this channel.
All in all we have 5N R + 20(+3) free parameters for γ i = 0 (γ i = 0). Clearly this number is larger than the number of parameters of two single-channel Breit-Wigner analyses, however, the advantage of the approach advocated here is that it allows for a combined analysis of all channels in a way that preserves unitarity, and for a straightforward inclusion of the 4π channel in the analysis. Note that the scalar resonances studied here are known to have prominent decays into four pions [5]; cf. also theoretical approaches modeling some of them as dynamically generated ρρ resonances [31][32][33][34].
The LHCb collaboration extracted two additional Swave resonances from their data [28], namely f 0 (1500) and f 0 (1790). Since there is no f 0 (1790) in the listings of the Review of Particle Physics by the PDG [5], we use the name f 0 (2020) for the higher state, in particular since the parameters we extract below are close to those reported for that resonance. The first fit includes our parametrization with N R = 2 and γ i = 0 (Fit 1). To test the stability of this solution, we also include a fit with N R = 2 and γ i = 0 (Fit 2) as well as a fit with N R = 3 and γ i = 0 (Fit 3). In order to obtain an estimate of the systematic uncertainty, we repeat each fit with two different assumptions about the third channel, which we take to be dominated by either σ σ or ρρ. The respective reduced χ 2 of the best fit results are listed in Table 1. We show the corresponding angular moments in Figs. 4 (ρρ) and 5 (σ σ ). In principle we could have also investigated mixtures of σ σ and ρρ intermediate states or parametrizations representing the channels π(1300)π or a 1 (1260)π reported to be relevant for the f 0 (1500) [5], however, since with the given choices we already find excellent fits to the data although the corresponding two-point function Σ 33 look vastly different for the σ σ and the ρρ case (cf. the lower right panel of Fig. 2), studying other possible decays will be postponed until data for further exclusive final states become available.
We note first of all that the ρρ fits have a lower reduced χ 2 compared to the σ σ fits. Allowing for a linear term in the source further improves the data description, as witnessed by the differences of Fits 1 and 2. The overall best reduced χ 2 is obtained by including another, third, resonance.
For the ρρ fit (see Fig. 4) we see that Fit 2 improves the description of Y 0 0 ππ in the energy region between 1.6 and 2.0 GeV. The biggest change between Fit 3 and the other ones is given by the better description of the high-energy tail in the decayB 0 s → J/ψ K + K − . For the σ σ fit, Fig. 5, we observe a similar picture. Fit 2 provides a very slight overall improvement of Fit 1. However, here the main difference between Fit 3 and the rest resides in the better description of the f 0 (1500) especially for the decaȳ B 0 s → J/ψ π + π − , while the high-energy tail ofB 0 s → J/ψ K + K − remains nearly untouched.
For a better comparison of the different fits we discuss the resulting form factors Γ s i in some detail. We begin by comparing the strange scalar pion form factor Γ s 1 as shown in Fig. 6. In all fits three resonances are clearly visible, namely the f 0 (980), f 0 (1500), and a broad structure around 2 GeV related to the f 0 (2020) resonance. Furthermore we also know that the input contains the broad f 0 (500) resonance. Fit 3 contains an additional resonance: in the case of the ρρ fit, it has its pole around 2.4 GeV and is relatively narrow. Notice that the maximum energy available for the ππ system in the decay studied is 2.27 GeV, thus this additional resonance in fact only contributes with its low-energy tail, giving small corrections for the high-energy parts of the angular moments. This is clearly visible in Y 0 0 K K at high energies in Fig. 4, where Fit 3 can describe the last data points better than Fits 1 and 2. In comparison we see that the σ σ fit lacks any such highenergy resonance. For this fit the difference between Fit 3 and the rest is only visible in the argument of Γ s 1 , showing a shift in the range 1.5 . . . 2 GeV. This improves the description of Y 0 2 ππ near the f 0 (1500) resonance. From this discussion it becomes clear that the data analyzed here do not allow us to extract information on any further resonance beyond f 0 (500), f 0 (980), f 0 (1500), and f 0 (2020).
By comparing the extracted kaon form factors Γ s 2 in Fig. 7 we see very similar features as for the pion form factor. However, the f 0 (1500) couples more weakly to the KK channel than to ππ, which is in line with what is reported about this state by the PDG [5]. The impact of the additional resonance in Fit 3 that appears outside the accessible data range is even more pronounced.
In Fig. 8 we compare the form factor of the additional, effective 4π , channel Γ s 3 . We see that the results of the fits with the 4π channel parametrized as ρρ differ significantly from the ones employing the σ σ variant. Moreover, also Fits 1-3 differ strongly from each other, even in the kinematic regime that can be reached inB 0 s decays. To further constrain these amplitudes it is compulsory to include data onB 0 s → J/ψ4π in the analysis, which is so far unavailable in partial-wave-decomposed form [44].
Finally in Fig. 9 we show the phases, δ, and inelasticities, η, that result for T 11 in the different fits, where we use the standard parametrization In the figure we also show the two-channel input phase δ 0 and inelasticity η 0 introduced in Eq. (10) as black solid lines. The comparison of the different lines demonstrates that the high-energy extension maps smoothly onto the low-energy input, as it should. In the phases one clearly sees the effect of the f 0 (1500), which leads to a deviation of the phase of T 11 from the input phase. In the inelasticity the full model starts to deviate from the input already at about 1.1 GeV as a consequence of the inclusion of the 4π channel. As in the phase the f 0 (1500) also leads to a pronounced structure in the inelasticity. It is interesting to observe that neither in the phase nor in the inelasticity there is a clear imprint of the f 0 (2020), which can be understood from its small coupling to the two-pion channel. In Fig. 9 we also show a comparison of our phases and inelasticities to those extracted in Ref. [45] (plotted as purple dashed lines) and the preferred solution [7] of the CERN-Munich ππ experiment [46] (data points with error bars). As one can see in the phase shifts, all analyses agree up to about 1.5 GeV. However, the effect of the f 0 (1500), present in all analyses, is very different. Also for the inelasticity there is no agreement between our solution and those from the two other sources, but here the deviation starts basically with the onset of theK K channel; for a more detailed discussion of the current understanding of the inelasticity in the scalar isoscalar channel, we refer to Ref. [11]. Note that there is also no agreement between the amplitudes of Ref. [7] and Ref. [45]. Thus, at this time one is to conclude that T 11 above 1.1 GeV is not yet known.
In a similar way, we can also compare the extracted ππ → KK amplitude T 12 with its absolute value g as well as its phase ψ, which are both shown in Fig. 10. While the resonance effects of the f 0 (1500) look qualitatively welldescribed by our high-energy extension, we see some differences to the actual data [47,48]. Note that the shown results are a prediction based solely on theB 0 s decay data and could be improved upon by explicitly taking the phase motion into account in the fit.

Extraction of resonance poles
In this section we present the extraction of resonance poles in the complex s-plane from the parametrizations discussed above. Traditionally those are given in terms of a mass M and a width Γ , connected to the pole position s p via [5] √ For narrow resonances located far from relevant thresholds, these parameters agree with the standard Breit-Wigner parameters. However, for broad and/or overlapping states, significant deviations can occur between the parameters derived from the pole location and those from Breit-Wigner fits. Since the analytic continuation to different Riemann sheets needs the on-shell scattering T -matrix as input, which, due to left-hand cuts induced by crossing symmetry, has a complicated analytic structure that cannot be deduced from the phase shifts straightforwardly, we use the framework of Padé approximants to search for the poles on the nearest unphysical sheets. For a thorough introduction into this topic, see e.g. Refs. [49][50][51].
As the form factor Γ s 1 (s) (Fig. 6) as well as T 11 (s) (Fig. 9) are smooth functions when moving from the upper complex s-plane of the first Riemann sheet to the lower complex s-plane of the neighboring unphysical sheet, we may expand both around some properly chosen expansion point s 0 according to The denominator allows for the inclusion of M resonance poles lying on the unphysical Riemann sheet. In the following we set M to 1, allowing for the extraction of the resonance that lies closest to the expansion point s 0 . The numerator ensures the convergence of the series to the form factor or the scattering matrix for N → ∞. In order to obtain the complex parameters a n and b n , we fit Padé approximants to both the form factor and the scattering matrix simultaneously. As both T 11 and Γ s 1 have the same poles, the parameters b n are the same for both, however, the a n are different. Note furthermore that the a 0 parameters are constrained by Γ s 1 (s 0 ) or T 11 (s 0 ), respectively.
For near-threshold poles such as the f 0 (500) and f 0 (980), we perform the Padé approximation not in s, but in the conformal variable instead [49]. This variable transformation maps the upper half complex s-plane of the first Riemann sheet to the inner upper half of a unit circle in the complex w plane, without introducing any unphysical discontinuities. The lower half of the second Riemann sheet is then mapped onto the lower half of the unit circle in the complex w-plane. This allows us to search for the two lowest poles within a circle around the expansion point s 0 , without being limited by the proximity of the ππ and KK thresholds, which are automatically taken care of. The statistical uncertainty is obtained through a bootstrap analysis of the fit results presented in Sect. 3.2. The systematic uncertainty coming from the Padé approximation on the other hand is estimated by [50] where s N p denotes the pole extracted by employing P N 1 (s, s 0 ). As in principle the results still depend on the expansion point s 0 , we proceed as follows. We first calculate Padé approximants for a varying s 0 ; near the true pole position, Corresponding residues of the poles are then described by the coupling strength g Rππ of the resonance R to ππ and the coupling g Rss of thess source to the resonance R. They are defined by the near-pole expansions [27,52] lim s→s p The extracted poles and residues for the resonances are shown in Table 2.
As we did not include any variation of the input phases, we see that the statistical uncertainty coming from the fit parameters of the higher-mass resonances has only a small impact on the poles of f 0 (500) and f 0 (980). In fact the uncertainty is dominated by the systematic error coming from the Padé expansion. At higher energies the statistical uncertainty becomes more significant.
However, overall we have strong systematic effects due to the assumptions on the parametrization such as the number of additional resonances and the linear terms in the polynomials. As we do not have a criterion that allows us to decide which fits we should prefer, we keep them all and perform a conservative estimate of the uncertainty: we choose a range for the resonance parameters such that all poles with their corresponding errors are included. The quoted mean is the middle of the resulting box as illustrated in Fig. 11.
In order to see whether the pole extraction leads to sensible results, we first compare our findings for the f 0 (500) and f 0 (980) to the literature [27,40,52,53]. In our parametrization the f 0 (500) has a mass of (442 ± 2) MeV with a width of (512 ± 10) MeV. For the f 0 (980) we find a mass of (996 ± 6) MeV and a width of (57 ± 11) MeV. As Ref. [40] serves as our input below 1 GeV, their pole positions are taken as a benchmark, which lie at (441 − i 544/2) MeV and (998 − i 42/2) MeV, respectively. While the real parts are therefore perfectly consistent, we see that our parametriza- Furthermore we can compare the coupling strengths g Rππ and g Rss to the ones found in Ref. [27], which we adjust for the fact that the latter are quoted for the complex conjugate poles. For the f 0 (500), we obtain g f 0 (500)π π = (4.53 ± 0.03) GeV, arg g f 0 (500)π π = (−73 ± 2) • , This is to be compared to g f 0 (500)π π = 4.76 GeV and arg g f 0 (500)π π = −76.4 • as well as g f 0 (500)ss = 17 ± 5 +1 −7 MeV and arg g f 0 (500)ss = 80.2 • [27]. With the exception of |g f 0 (500)π π |, which appears to be shifted by about 5%, these numbers are consistent with our findings. For the f 0 (980) pole, we find g f 0 (980)π π = (3.1 ± 0.5) GeV, arg g f 0 (980)π π = (−81 ± 5) • , in comparison to the reference values g f 0 (980)π π = 2.80 GeV, arg g f 0 (980)π π = −85.3 • , g f 0 (980)ss = 146 ± 44 +14 −7 MeV, and arg g f 0 (980)ss = 14.2 • [27]. In this case therefore all parameters are consistent within uncertainties, with a small tension for the argument of g f 0 (980)ss . In particular, we reproduce the well-known hierarchy in the couplings to thess current: the f 0 (980) couples to the strange scalar current an order of magnitude more strongly than the f 0 (500) does. Overall we find good agreement of our pole parameters for f 0 (500) and f 0 (980) with the literature. We see that, a posteriori, the subtraction of the additional term in the scattering amplitude that introduces the explicit resonances, cf. Eq. (22), suppresses its influence on the  [45]. In addition we plot the preferred phase shifts and inelasticities [7] of the CERN-Munich ππ experiment [46], which are denoted by data points with error bars lower-mass poles sufficiently. The agreement between the reference parameters and ours gives us confidence for an extraction of the higher poles via Padé approximants. As a reference for the higher resonance poles, we compare to the Breit-Wigner parameters of LHCb [36]. For the f 0 (1500), the collaboration quotes a resonance with mass (1465.9 ± 3.1) MeV and width (115 ± 7) MeV. The pole we extract corresponds to a mass of (1465 ± 18) MeV and a width of (100 ± 19) MeV, which lies within the previously quoted uncertainties of LHCb. The uncertainties we find are significantly larger: this is most likely due to the more flexible range of resonance parametrizations we employ; the masses and widths extracted using Breit-Wigner functions only are probably too optimistic. In addition we can extract the corresponding residues, which are given by g f 0 (1500)π π = (2.9 ± 1.0) GeV, arg g f 0 (1500)π π = (−42 ± 4) • , g f 0 (1500)ss = (125 ± 76) MeV, arg g f 0 (1500)ss = (167 ± 21) • .
The main uncertainties stem from the assumptions made on the parametrization of the form factor, such as the number of resonances and the additional channels. Nevertheless, we note that, despite a large uncertainty, the central value for g f 0 (1500)ss seems to be comparable to g f 0 (980)ss . For further comparison, according to Refs. [27,54] the a 0 (1450) couples to an isovector scalarūd current with g a 0 (1450)ud = (284 ± 54) MeV, which is of the same order as our extracted value for g f 0 (1500)ss . The precise relation between the two couplings might be used to elucidate the structure of a scalar nonet around 1.5 GeV, which is however beyond the scope of the present study. For broad, overlapping resonances a definition of branching ratios is not straightforward. Here we follow a prescription originally proposed to define the width of f 0 (500) → γ γ [55] by using the narrow-width formula of the form is shown in blue, Fit 2 in red, Fit 3 in green, and the input g 0 and ψ 0 [40] in black. For comparison we show the amplitude analyses of Refs. [47] (open diamonds) and [48] (filled stars) with the residues as coupling constants. With this we can deduce a branching ratio B f 0 (1500)→ππ = (58±31)%, where the main uncertainty stems from the difference between Fits 1 and 2 with an additional σ σ channel compared to the rest of the fits. This is compatible with the (much more precise) branching ratio quoted by the PDG, B f 0 (1500)→ππ = (34.9 ± 2.3) % [5]. The last resonance identified by LHCb as the f 0 (1790) has a mass of (1809 ± 22) MeV with a width of (263 ± 30) MeV. As we do not impose a Breit-Wigner line shape, our fits seem to prefer a significantly heavier and much broader resonance with mass (1910 ± 50) MeV and a width of (398 ± 79) MeV. Note that for the average we neglected the pole extracted from Fit 1 with the ρρ parametrization, since this fit describes the prominent resonance structure in the ππ spectrum less accurately than the rest of the fits. As the pole position of the higher pole extracted in our analysis is in better agreement with the f 0 (2020) of the PDG (which quotes a mass of (1992 ± 16) MeV and a width of (442 ± 60) MeV [5]), we will refer to it as such in the following. Furthermore we see that this pole allows for a stronger variance in the different fits. As its line shape does not only depend on the interference with other resonances, but also on further inelasticities, additional information about these channels would be appreciable.
As we can see the coupling strength to the ππ-channel is consistent with 0 within 1.5σ . The big uncertainty also strongly influences the extraction of g f 0 (2020)ss , which in addition is affected by a strong systematic uncertainty coming from the parametrization and can hardly be constrained in a meaningful manner. Using the narrow-width formula of Eq. (44), the branching ratio into ππ is B f 0 (2020)→ππ = (1.3 ± 1.8)%, which is obviously also consistent with zero. No meaningful branching ratios are quoted by the PDG in this case.
Since the bare resonance coupling strengths g r i as well as the bare resonance masses m r are source-independent, Table 2 Padé poles for f 0 (500), f 0 (980), and f 0 (1500) for N = 5, as well as f 0 (2020) for N = 6. The error is the uncorrelated sum of statistical and systematic uncertainty