Insights into $Z_b(10610)$ and $Z_b(10650)$ from dipion transitions from $\Upsilon(10860)$

The dipion transitions $\Upsilon(10860)\to\pi^+\pi^-\Upsilon(nS)$ ($n=1,2,3$) are studied in the framework of a unitary and analytic coupled-channel formalism previously developed for analysing experimental data on the bottomoniumlike states $Z_b(10610)$ and $Z_b(10650)$ [Phys. Rev. D 98, 074023 (2018)] and predicting the properties of their spin partners [Phys. Rev. D 99, 094013 (2019)]. In this work we use a relatively simple but realistic version of this approach, where the scattering and production amplitudes are constructed employing only short-ranged interactions between the open- and hidden-flavour channels consistent with the constraints from heavy quark spin symmetry, for an extended analysis of the experimental line shapes. In particular, the transitions from the $\Upsilon(10860)$ to the final states $\pi \pi h_b(mP)$ ($m=1,2$) and $\pi B^{(*)}\bar B^* $ already studied before, are now augmented by the $\Upsilon(10860)\to\pi^+\pi^-\Upsilon(nS)$ final states ($n=1,2,3$). This is achieved by employing dispersion theory to account for the final state interaction of the $\pi\pi$ subsystem including its coupling to the $K\bar K$ channel. Fits to the two-dimensional Dalitz plots for the $\pi^+\pi^-\Upsilon$ final states were performed. Two real subtraction constants are adjusted to achieve the best description of the Dalitz plot for each $\Upsilon(nS)$ $(n=1,2,3)$ while all the parameters related to the properties of the $Z_b$'s are kept fixed from the previous study. A good overall description of the data for all $\Upsilon(10860)\to\pi^+\pi^-\Upsilon(nS)$ channels achieved in this work provides additional strong support for the molecular interpretation of the $Z_b$ states.


I. INTRODUCTION
The spectroscopy of hadronic states containing heavy quarks remains one of the fastest developing and most intriguing branches of strong interaction studies. Many new states have been discovered in the spectrum of charmonium and bottomonium which do not fit into the quark model scheme and qualify as exotic states. For example, the states Z ± b (10610), Z ± b (10650) [1], Z ± c (3900) [2,3], Z ± c (4020) [4], Z ± (4430) [5][6][7][8] are charged and decay into final states containing a heavy quark and its antiquark. Since the production of this pair of heavy quarks in the decay is highly suppressed in QCD, it must have been present in the wave functions of the states and as such the Z Q states cannot be conventionalQQ (with Q denoting a heavy quark) mesons as their minimal quark contents is four-quark. The interested reader can find a comprehensive overview of the current experimental and theoretical status of the exotic hadrons with heavy quarks in dedicated review papers, for example, in Refs. [9][10][11][12][13][14][15].
The Z ± b (10610) and Z ± b (10650) bottomoniumlike states (in what follows often referred to as Z b and Z b , respectively) are ideally suited for both experimental and theoretical studies since there exist two resonances in the same J P C = 1 +− channel [16] split by only about 40 MeV which are simultaneously seen in several modes. Specifically, the Belle Collaboration observed them as distinct peaks (i) in the invariant mass distributions of the π ± Υ(nS) (n = 1, 2, 3) and π ± h b (mP ) (m = 1, 2) subsystems in dipion transitions from the vector bottomonium Υ(10860) [1] and (ii) in the elastic BB * 1 and B * B * channels in the decays Υ(10860) → πB ( * )B * [17,18]. The two most prominent explanations for the Z b 's claimed to be consistent with the data are provided by a tetraquark model [19][20][21] and a hadronic molecule picture [22][23][24][25][26][27][28][29][30][31]. A review of the sum rules approach to the exotic states with heavy quarks and relevant references on the subject can be found in a recent review [32]. It should be noted that a particularly close location of the Z b 's to the thresholds of the BB * and B * B * channels, which in addition are the most dominant decay modes for them, provides a strong hint in favour of their molecular interpretation.
Both the Z ± b (10610) and Z ± b (10650) contain a heavy bb pair, so it is commonly accepted that the heavy-quark spin symmetry (HQSS) should be realised to high accuracy in these systems and indeed, HQSS is able to explain naturally the interference pattern in the inelastic channels Z ( ) b → πΥ(nS) and Z ( ) b → πh b (mP ) [22]. In Ref. [33], an effective field theory (EFT) approach to the Z b states consistent with HQSS and chiral symmetry was developed to perform a combined analysis of the experi-mental data in the channels The information on the branching fractions in the transitions Υ(10860) → πZ ( ) b → ππΥ(nS) (n = 1, 2, 3) was also used, but no analysis of the line shapes in these channels was performed for the reasons explained below. A fairly good description of the data was achieved in different fitting schemes described in detail in Ref. [33]. As expected, the experimental data on the Z b 's are fully consistent with HQSS, since symmetry violating terms in the effective hadronic potential are argued to play a minor role [33]. In Ref. [34], the approach was extended to predict in a parameter-free way the properties of the spin partner states of the Z b 's, the W bJ 's (J = 0, 1, 2).
The one-pion exchange (OPE) in the bottomoniumlike systems under consideration was a special concern of the quoted works [33,34], and it was concluded to play an important role for the Z b 's and W bJ 's. Indeed, the poles of the Z b 's and W bJ 's that were originally classified as virtual states in the pionless framework moved above the nearby elastic thresholds to become resonances, as an effect of the OPE. Meanwhile, the conclusion that all these states are hadronic molecules, based on a decent description of the data, follows already from the scheme with purely contact interactions in the B ( * )B * system (Scheme A, in the notation of Ref. [33], yields χ 2 /N dof ≈ 1.23). Note also that this fitting scheme provides results identical to those obtained with the help of an analytical parametrisation for the line shapes derived previously in Refs. [35,36].
Not all experimental information used in the aforementioned combined analysis could be considered on equal footing. Indeed, while the line shapes in the πh b (mP ) and B ( * )B * channels could be fitted directly, as discussed above, only the total branchings for the πΥ final states were used in the fit. The signal in the latter channels contains a significant nonresonant contribution that depends on the invariant mass of the two-pion system, so that the amplitude analysis has to be multidimensional. This analysis is in the spotlight of the present work. In particular, we generalise the approach developed in Ref. [33] to incorporate coupled-channel effects from the ππ-KK interactions in the final state using a model-independent dispersive approach. Then we perform maximum likelihood fits to the Dalitz plots of the reactions Υ(10860) → ππΥ(nS) (n = 1, 2, 3). To keep consistency with the data in the πh b (mP ) and B ( * )B * channels, we directly employ the inelastic production amplitudes obtained in Ref. [33] for Scheme A as input for the present research. Since the focus of the present study is on the development of the dispersive treatment of the final state interactions (FSI), we resort to a simple pionless formulation, as provided by Scheme A, while effects from the OPE will be included in future studies. Thus in this study we focus on the following goals: (a) A development of a dispersive approach to the Υ(10860) → ππΥ(nS) transitions and a systematic account for the effects from the ππ FSI including the coupling to the KK channel. While for the Υ(2S) and, especially, Υ(3S) in the final state the ππ-KK coupling is expected to play a marginal role, it should be important for the Υ(1S) channel near the KK threshold. This effect can be included in a model-independent way using an Omnès matrix constructed from high accuracy determinations of the ππ and KK scattering amplitudes as well as from the B s decay data [37,38].
(b) Our focus is on the inclusion of the FSI while keeping the full complexity of the Z b dynamics, so we consider two production mechanisms for the transitions Υ(10860) → ππΥ(nS), namely (i) through the contact operators with two real parameters and (ii) through B-meson production assuming pointlike vertices with the subsequent B-meson interactions in the final state, that is, via the process Υ(10860) → B ( * )B * π → ππΥ(nS). Both mechanisms are supplemented with the ππ FSI. Note that in Ref. [39] also a possible impact of the boxdiagram mechanism was studied, which is not included here. The underlying rationale is that in the Υ(10860) decays the Z b states can go on-shell and should by far dominate the effects from the B ( * ) B * intermediate states. The corresponding imaginary parts in the production amplitudes are taken into account explicitly in this work. As a consequence, only two real subtraction constants defined in the mechanism (i) are sufficient to dispersively reconstruct the amplitude, which is insensitive to the high-energy integration range. This is unlike to Ref. [41], where two complex coefficients were utilized in a related study of the dipion transitions in the charmonium sector.
(c) The Dalitz plots for the Υ(10860) → ππΥ(nS) transitions contain nontrivial information about the Z b 's -these states can be clearly seen in the πΥ(nS) invariant mass distributions and have imprint also on the ππ spectrum. Thus we analyse the two-dimensional Dalitz plots to check whether the results for the Z b 's from our previous analyses are consistent with them.
The paper is organised as follows. In Sect. II we briefly introduce the coupled-channel approach suggested and used in Refs. [33,34]. In Sect. III a dispersive approach to the decay amplitude is developed to take into account the ππ interaction in the final state. Section IV is devoted to the data analysis for the reactions Υ(10860) → ππΥ(nS) (n = 1, 2, 3). Our conclusions are discussed in Sect. V. Appendices A and B provide some technical details of the dispersive approach used in this work, including a discussion of the anomalous contributions to the amplitude.

II. COUPLED-CHANNEL APPROACH
In this section we briefly recall some essentials of the coupled-channel approach previously developed in Ref. [33] to perform a combined analysis of the data for the bottomoniumlike states Z b (10610) and Z b (10650). The channels with hidden bottom (labelled by latin letters), (2) are referred to as inelastic ones while the open-bottom channels (labelled by greek letters), are denoted as elastic ones. The interaction potential between different channels takes the form of a matrix, The main purpose of the present work is to incorporate into the current coupled-channel scheme of Ref. [33] the pion interaction in the ππΥ(nS) final states. Consequently, although the OPE was argued in Ref. [33] to provide an important contribution to the elastic potential, its analytic structure is quite complicated and had the largest impact on the BB * channel. Accordingly, it will be neglected in the current study, especially since existing data can be quite well described within the purely contact Scheme A of Ref. [33]. Thus, in what follows, we stick to this scheme and assume that the leading left-hand cut contributions to the amplitude for the ππ FSI are generated by the Z b 's poles and the B ( * )B * cuts which are taken into account in the present approach. To be specific, in this work we employ the following approximations: • Only O(p 0 ) contact interactions are included in the elastic channels, so that the v αβ (p, p ) from Eq. (4) takes the form where C d and C f are independent low-energy constants; • Elastic-to-inelastic transition potentials are parametrised via coupling constants as where k i and l i are the momentum and the angular momentum in the i-th inelastic channel, respectively. The inelastic momentum is calculated as where m Hi (m hi ) is the mass of the heavy(light) meson in this channel, M is the total energy of the system, and λ(m 2 1 , m 2 2 , m 2 3 ) is the standard Källen triangle function, λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2xz − 2yz. (8) The coupling constants g iα are constrained by HQSS: = −1, n = 1, 2, 3, Therefore, as in Ref. [33], they will be quoted only for the BB * channel in the form • Following the arguments from Refs. [35,36], direct interactions in the inelastic channels are neglected, v ji (k j , k i ) = 0.
• The long-ranged part of the pion exchange between the B ( * ) mesons is not considered 2 .
As a result of these approximations, the effective elastic-to-elastic channel transition potential takes the form where U α (M, p) denotes the physical production amplitude of the α-th elastic channel from a pointlike S-wave  [33]. The cut-off Λ is set to 1 GeV (see the discussion in the quoted paper). For the inelastic coupling constants only the absolute values are presented since physical quantities are not sensitive to their signs.
source, and F BB * (M, p) = −F B * B * (M, p) = 1 as dictated by HQSS. The Green's function for a two-heavymeson intermediate state reads where m α th stands for the α-th elastic threshold and µ α is the reduced mass in this channel. Other components of the multichannel amplitude responsible for production of the inelastic channels in the final state can be obtained from U α (M, p) algebraically, which is a consequence of the omitted direct interactions in the inelastic channels. In particular, for the i-th inelastic channel in the final state we have with the momentum k i defined in Eq. (7) above. It has to be noticed that the Born amplitudes F i (M, p) coming from the inelastic sources were neglected in Eq. (14). This is justified for the πh b (mP ) channels, where the data are dominated by the Z b (10610) and Z b (10650) poles emerging from the B ( * )B * dynamics. The corresponding line shapes were included into the combined fit performed in Ref. [33]. On the contrary, in the heavy-spinconserving πΥ(nS) channels, the Born term needs to be kept and the ππ interaction in the final state has to be included. How this can be done in a model-independent way will be discussed in detail below.
The one-dimensional distributions for the differential widths in the elastic (B ( * )B * ) and inelastic (πh b (mP )) channels used in Ref. [33] read respectively, where p * π is the three-momentum of the spectator pion in the rest frame of the πΥ(10860) and p α (p i ) is the three-momentum in the α-th elastic (ith inelastic) channel in the rest frame of the B * B( * ) (πΥ(nS)/πh b (mP )) system. Then, the total branching fraction in an elastic or inelastic channel x is defined as where and the integral in M covers the entire kinematically allowed region for the considered channel x.
The line shapes obtained in the fitting scheme described above are presented in Fig. 1, and the parameters of the fit are listed in Table I. As was explained above, only the total branchings of the πΥ channels were included in the fits.
In what follows, the left-hand cut structure of the multichannel production amplitude U i (i = πΥ(nS)) from Eq. (14), obtained in the framework of the contact Scheme A, will be used as input for a dispersive reconstruction of the ππ FSI.

A. Kinematics of the reaction
In this subsection we introduce the kinematics of the decay Υ(10860)(p i ) → Υ(nS)(p f )π + (p 1 )π − (p 2 ) with n = 1, 2, 3. Following a standard approach to such reactions, we built the amplitude M (s, t, u) in a crossed channel, Υ(p i ) + Υ (p f ) → π(p 1 ) + π(p 2 ), and define the Mandelstam invariants accordingly, where m π , m i , and m f are the masses of the pion, Υ(10860) ≡ Υ, and Υ(nS) ≡ Υ , respectively. Thus, In order to proceed, we resort to the kinematics in the centre-of-mass frame of the two pions in the final state, so that (z ≡ cos θ, where θ is the angle between the 3momenta p 1 and p f ), with the function λ defined in Eq. (8) above. Consequently, z can be expressed in terms of t and u as Since the production amplitude for the process Υ(10860) → π + π − Υ(nS) has the form 3 the double differential production rate can be written as where the overall normalisation constant N will be fitted to the data.
In this subsection we introduce the meson-meson interaction in the final state. The partial wave decomposition of the amplitude M (s, t, u) reads where P l (z) are the Legendre polynomials and the sum runs over all relevant angular momenta l. We start from the amplitude M (s, t, u) projected onto the ππ S-wave, which can be split into two pieces, where the first and second term contain the right-and left-hand cuts only, respectively. The right-hand cut of the amplitude M R 0 comes from the FSI while the lefthand cuts of the amplitude M L 0 are due to the dynamics related to the Z b states. If the contribution M L 0 is known and only the ππ channel is considered for the FSI, the full amplitude can be reconstructed dispersively via the solution of the inhomogeneous Omnès problem as (see Ref. [42] for a related discussion) where Ω 0 (s) is the S-wave single-channel Omnès function 4 and δ is the ππ S-wave phase shift (see Appendix A for details). However, given that the energy in the ππ system in the reaction Υ(10860) → ππΥ(1S) extends to 1.4 GeV, that is far beyond the KK threshold, the inclusion of the KK component becomes necessary. Generalisation of Eq. (28) to multiple channels is straightforward, Here, the multichannel Omnès matrix obeys the matrix equationΩ 4 We use the standard notation Ω I l for the Omnès function where l and I stand for the partial wave and isospin, respectively. However, since in this work we deal only with isoscalars, the superscript I = 0 is omitted everywhere.
where hats indicate multicomponent objects (vectors and matrices),σ(s) = diag{σ π , σ K } is a diagonal matrix with σ P (s) = 1 − s th P /s, and s th P for the threshold in the corresponding channel (P = π, K). In particular, we have wave meson-meson coupled-channel amplitudeT can be parametrised by the the ππ scattering phase shift δ(s) [43][44][45][46] as well as the absolute value and phase of the ππ → KK transition [45,46], g(s) and ψ(s), respectively, as To get the two-pion FSI amplitude one has to consider the component [M 0 (s)] ππ of the vector (29). If the amplitude contains contributions from higher partial waves while the FSI is taken into account only in the S wave, one can writê whereM no-FSI =M L 0 +M higher is the complete tree level production amplitude in the t-and u-channel, not projected onto partial waves, while the effect of the FSI is taken into account by the second term in Eq. (32). In this study, the ππ component of the production ampli-tudeM no-FSI and its S-wave projectionM L 0 are adopted from Ref. [33] -see Sec. III C for a detailed discussion. Meanwhile, the resonance production in the channel Υ(10860) → KKΥ(nS) which proceeds through the B-and B s -meson loops is not considered since no information about the SU (3) partners of the Z b states is available yet.
The dispersive integral in Eq. (32), where the lower index indicates l = 0 for the S wave, may need to be subtracted n times to improve convergence and to diminish the role played by the large-s region where the ππ scattering phase is not known well enough. Then, one arrives at whereP n−1 (s) is a polynomial of the order n − 1. If the amplitudeM L 0 (s) has both real and imaginary parts, then the polynomial coefficients are complex numbers. Meanwhile, if there are good reasons to believe that the imaginary part of the amplitude ImM L 0 (s) is controlled by well understood physics (see also a related discussion in Sec. III D below), then the imaginary part of the polynomial P n−1 (s) can be evaluated exploiting sum rules via where it was used that the quantityΩ −1 0 (s )T (s )σ(s ) is real. This allows one to re-write Eq. (34) in the form where the polynomial R n−1 (s) is real by construction, and so are its coefficients. As discussed below, ImM L 0 (s) is non-vanishing on a finite interval of s only and, accordingly, the integral in the last line of Eq. (36) does not require any subtractions.

C. The left-hand cut production amplitude
In order to proceed with the formulae derived in the previous subsection, we need to specify the form of the production amplitude M no-FSI introduced in Eq. (32) and determine its S-wave projection, Consider, as a preliminary step, a stable-Z b exchange in the t-and u-channel. Then, assuming pointlike Υ ( ) → πZ b vertices, up to an overall constant, one can write the invariant Born amplitude as where m z is the mass of the mentioned stable Z b particle. For future convenience we specify this mass as an argument of the amplitude.
Performing the partial wave projection as introduced in Eq. (38) for the Born amplitude (39) Here is the root of the equation Y (s) = 0. Furthermore, the logarithmic branch points s ± (also known as anomalous thresholds) found as the roots of the equation t(s ± , z = ±1) = m 2 z read where λ is the triangle function from Eq. (8).
In the regime s + is real and the anomalous threshold generates only a phase term which is included in the first formula in Eq. (40) (see also Ref. [41] for a related discussion). However, for the branch point s + becomes complex and the dispersive integral defined in Eq. (36) acquires an additional anomalous contribution calling for an integration along some complex path (see Appendix B for details). Namely, using M L 0;anom = −4πi/κ for the anomalous discontinuity, the integralÎ and ζ = (1 − x)s + + x 4m 2 π is the straight-line path between the two-pion threshold and the branch point of the logarithm, s + .
A crucial point of the coupled-channel approach developed in Ref. [33] is that the resonances Z b are not introduced as asymptotic states of the theory but appear as near-threshold poles of the amplitude fitted to the data. This implies that, instead of the stable Z b propagator used in Eq. (39), the inelastic amplitude U i (i = πΥ(nS) with n = 1, 2, 3) from Ref [33], generated through the Bmeson loops and evaluated as given in Eq. (14), provides the input for building M no-FSI and M L 0 -see Fig. 2 for its diagrammatic representation. To proceed, we employ a dispersive representation for the production amplitude U (to simplify notations we omit the inelastic index i and thus consider a particular inelastic final state), where we used Eq. (39) and introduced the spectral func-tion ρ(µ 2 ) = − 1 π Im U (µ 2 ).
The lower limit in the integral above is given by the lowest relevant threshold which may contribute to the imaginary part of the amplitude. The unitarity of the production amplitudes U (µ 2 ) requires integration from the lowest inelastic threshold πΥ(1S), although the leading contributions start from the BB * π threshold. The limit of a stable particle, Eq. (39), is reached from Eq. (47) for ρ(µ 2 ) = δ(µ 2 − m 2 z ). With the mass distribution of the Z b states included, the S-wave partial wave amplitude reads where M L 0,stable (s, µ) is defined in Eq. (40). Also, the anomalous term from Eq. (46) has to be weighted with the spectral function to read

D. Matching to chiral perturbation theory
Two comments on the convergence of the subtracted dispersive integral (36) entering Eq. (32) are in order here. First, as follows from Eqs. (40) and (43), the imaginary part of the M L 0 (s) may be nonzero only at a finite interval in s between the logarithmic branch points s − and s + . Therefore, the integral from ImM L 0 in the last term in Eq. (36) is finite. The appearance of an imaginary part in the left-hand cut amplitudes is a specific consequence of the cuts which are allowed in the process Υ(10860) → ππΥ(nS), especially from the B ( * )B * π intermediate states but also from inelastic channels. While inelastic cuts can also contribute to similar decays from the Υ(3S) and Υ(4S), the presence of the cuts in the elastic channels is only allowed kinematically starting from the Υ(10860).
The number of subtractions n needed to render the dispersion integral convergent, can be determined by the high-energy behaviour of the function Ω −1 0 (s )T (s )σ(s )ReM L 0 (s ). Since the Omnès matrix is determined by an unsubtracted dispersion relation (30), each of its elements behaves as 1/s at high energiessee Ref. [47] for details. Thus each element of its inverse scales as s. For the scattering amplitudeT in the twochannel case defined in Eq. (31) it is possible to demonstrate that for T 12 ∝ 1/s 3/2 at large s [48], T 11 and T 22 need to scale as 1/s 3 . Therefore, we proceed with a conservative estimate thatT (s) ∝ 1/s 3/2 at large s. The lefthand cut amplitude for a stable particle M L 0,stable (s, m z ) falls off as log(s)/s at large s, however the full left-hand cut production amplitudeM L 0 may decrease slower and is expected to approach a constant. Then, even without subtractions, one in principle should arrive at a convergent integral. However, to suppress the contribution of the large-s region, where the details of the ππ interaction are badly known, in what follows, twice subtracted dispersive integrals are considered, and the polynomial in Eq. (36) takes the form with real parameters a and b, as was explained above. In what follows, it will be shown that one of these constants is mostly redundant at least for production of Υ(nS) with n = 2 and 3.
It is important to notice that the polynomial R 1 (s) parametrises the amplitude for ΥΥ ππ at small values of s and as such can be matched to chiral perturbation theory. Specifically, in the limit of switching off the finalstate interactions, δ(s) → 0, g(s) → 0 in Eq. (31) and thus setting Ω 0 (s) → 1, the subtraction functions must agree with the chiral amplitudes corresponding to the direct transitions Υ → ππΥ [49].
If one introduces spin multiplets for heavy-heavy fields, then the effective Lagrangian for the contact ΥΥ ππ and ΥΥ KK coupling, at the lowest order in the chiral and heavy-quark expansions, reads [39,49,50] where v µ is the 4-velocity of the heavy quark. The contribution of the pseudoscalar Goldstone bosons for the spontaneous breaking of the chiral symmetry can be parametrised as  where q is the 3-momentum of the final Υ in the rest frame of the initial Υ, that is, Up to some small corrections, the amplitude (54) behaves as a linear polynomial in s. Thus the chiral amplitude at low energies depends on the two low-energy constants (LECs) c 1 and c 2 which can be treated as fitting parameters instead of a and b from Eq. (51). This amplitude corresponds to the contact diagram depicted in Fig. 3(a). Then, the amplitude M (s, t, u) from Eq. (24), which now includes the effects from the ππ and KK where M L s (s ) andÎ anom 0 (s) are given in Eqs. (47) and (50), respectively.
The diagrams representing different contributions to the amplitude of Eq. (56) are depicted in Fig. 3: the sum of the diagrams (a) and (b) corresponds to the termΩ 0 (s)M χ,ππ 0 (s), the diagram (c) gives the production amplitude M no-FSI in the t-and u-channel, and the diagram (d) describes the contribution of the last term Ω 0 (s)Î  where P 2 (z) is the second-order Legendre polynomial (see also Eq. (23)), the amplitude M χ,ππ 2 (s) extracted from the Lagrangian (52) reads and the diagrams which correspond to the amplitude (60) coincide with those depicted in Fig. 3 (a) and (b), however with no kaons in the loop. No additional parameter is involved in the amplitude (60), since c 2 also enters Eq. (54). The D-wave Omnès function Ω 2 (s) in Eq. (60) is calculated using the D-wave ππ phase shift from Ref. [43] and is dominated by the f 2 (1270) resonance contribution. In general, the amplitude in Eq. (59) should also contain the dispersive integralĨ (n2) 2 (s), which is however neglected in the current study. This is motivated by the fact that the corresponding D-wave contribution from the chiral polynomial in Eq.(60) plays only a very minor role in the fits, as discussed in Sec. IV. While in this study the main focus is put on the development of the appropriate formalism and testing the general consistency of the coupled-channel EFT approach of Ref. [33] with the data in the Υ(10860) → ππΥ(nS) decays, we (s) to a future global analysis of all data available in various channels.

IV. DATA ANALYSIS
With the coupled-channel approach developed in Ref. [33] and further augmented by the ππ/KK interaction in the final state, as explained in the previous section, we are in a position to analyse the data on the decays Υ(10860) → ππΥ(nS) (n = 1, 2, 3). We use the Belle data from Ref. [52] and make a maximum likelihood fit to the two-dimensional distributions. We exclude the region near the Dalitz plot boundary to minimise the effects of the e + e − centre-of-mass energy spread and the detector resolution. The signal function is corrected for the efficiency, and the experimental background distribution is added. To account for the invariant mass resolution in M 2 (πΥ(nS)), we convolute the theoretical results with the Gaussian resolution function with σ = 4. The results of the data analysis performed in this work are presented in Fig. 4 in the form of one-dimensional projections of the Dalitz plots, namely the invariant mass distributions M 2 (ππ) and M 2 (πΥ(nS)). The parameters of the fits are listed in Table II, where the last column shows that at least for production of Υ(2S) and Υ(3S) only a linear combination of the two constants is relevant. From Fig. 4 one can see that the developed approach is able to describe the data in the ππΥ channels rather well. In each plot we give several curves showing different contributions to the total rate. The relative importance of these contributions changes significantly with the mass of the bottomonium Υ(nS) in the final state.
• In case of the ππΥ(3S) channel, given a very limited phase space available, the distributions are dominated by the t-and u-channel production amplitudes from Ref. [33] (M no-FSI in Eq. (32)) which capture the gross features of the experimental distributions. Effects from the ππ FSI and chiral polynomials are marginal individually and, in addition, numerically cancel each other to a large extent.
• For the ππΥ(2S) channel, the resulting contribution from the t-and u-channel amplitudes and the dispersively reconstructed s-channel FSI are all important but not sufficient to explain the data. The missing strength and the energy dependence comes from the chiral contact terms with the ππ FSI.
• Finally, in case of the ππΥ(1S) final state with the maximal available phase space, the pattern of the individual contributions is qualitatively similar to that in the Υ(2S) channel. In particular, the t-and u-channel production amplitudes provide the peaking structures from the Z b 's and enhance the small ππ region in the corresponding line shape. The dispersive integral I 0 , enhanced by the ππ-KK FSI from the coupled-channel Omnès function, is very important. It provides more than half of the full signal at small M 2 πΥ(nS) and large M 2 ππ but also significantly affects the M 2 πΥ(nS) line shape in the Z b 's area. In addition, I 0 drives the shape of the M 2 ππ spectrum near the KK threshold and together with the t-and u-channel production amplitudes describes the low M 2 ππ region. The residual contribution comes from chiral contact terms in S wave, while their effect in D wave is minor.
In Fig. 5 we give the helicity angular distributions, including individual contributions to it. The anisotropies of these distributions are largely driven by the higher partial-wave contributions from the amplitude M no-FSI (t, u) in Eq. (59) -as usual the partial wave expansion in some subsystem converges badly, as soon as there is a narrow, near on-shell state in the crossed channel.
While the approach advocated in this study allows for a quite reasonable quantitative understanding of the line shapes in the ππΥ(nS) channels, we would like to emphasise that the peaks of the Z b 's in our results (red solid curves in Fig. 4), which by construction are consistent with the experimental B ( * )B * and πh b (mP ) (m = 1, 2) distributions, are not exactly in accord with the data in the πΥ(nS) channels. This observation still calls for an explanation.

V. CONCLUSIONS
In this work we extended the coupled-channel approach developed previously in Ref. [33] to incorporate the ππ FSI in the heavy-quark-conserving channels ππΥ(nS) (n = 1, 2, 3). Maximum likelihood fits to the twodimensional Dalitz plots are performed in each such channel. Compatibility with the previous analysis of the line shapes in the B ( * )B * and πh b (mP ) (m = 1, 2) channels is guaranteed by employing the amplitudes obtained in Ref. [33] as input for the current research.
The ππ FSI was incorporated into the coupled-channel scheme employing a dispersive approach, in which the left-hand cut contributions were provided by the previ-ously found inelastic production amplitudes. Also, it was conjectured that the dominating sources of the imaginary parts of the production amplitudes in the ππΥ(nS) channels are fully under control and, therefore, only the real parts of the complex polynomials were fitted to the data, while their imaginary parts appeared as actual predictions of the approach. As a cross-check, it was verified that the fits did not improve much if the imaginary parts of the chiral polynomials were also fitted to the data.
The results obtained demonstrate that the extended approach developed in this work is able to describe the existing experimental data with a reasonable accuracy. As expected, the role of the FSI increases with the increase of the allowed phase space. For example, it is remarkable, that the results for the ππΥ(3S) channel, which are almost completely driven by the production operators from Ref. [33] with all the parameters fixed from other data, are in a relatively good agreement with the data for the ππ and πΥ(3S) projections. Meanwhile, effects from the ππ FSI in the S-wave are found to be important for the ππΥ(2S) and, especially, ππΥ(1S) channel. In the latter case, also the KK channel plays a very important role. In particular, the coupled-channel ππ-KK dynamics shows up as a very clear dip structure near the KK threshold, which is consistent with the data. The effect of the D-wave resonance f 2 (1270) included via the D-wave Omnès function together with the contact production operators is found to be very minor for all channels.
Among the open questions to mention is an observation that the peaks corresponding to the Z b 's from the presented analysis do not completely correspond to the data in the πΥ(nS) channels. To address this question a full combined fit to all measured production and decay channels for the Z b 's should be performed. This lies beyond the scope of the present research.
We consider the results presented here as an important step towards a comprehensive combined analysis of the data in all production and decay channels of the Z b 's, which should in the future also incorporate pion exchanges. Such an analysis which will include the full information contained in the two-dimensional distributions for the ππΥ(nS) channels should provide the most accurate determination of the parameters of the theory and, as a result, allow making quite precise predictions for the line shapes, pole positions, decays couplings and other properties of the yet unobserved spin partner states W bJ . Such an insight is expected to provide an important theoretical background for the searches of the W bJ 's as well as other near-threshold exotic candidates in the experiment Belle II and, possibly, in LHCb and future hadronic experiments. The solution of Eq. (A4) for M R 0 can be found in the form where G(s) is an unknown function to be restored dispersively from its discontinuity. To this end we find from Eq. (A5) that where X ± ≡ X(s ± i ) for an arbitrary function X(s), and it was used that Ω 0+ ≡ Ω 0 = |Ω 0 |e iδ , Ω 0− = Ω 0+ e −2iδ = |Ω 0 |e −iδ . If the amplitude M L 0 (s) has no imaginary part, it is easy to verify that the phase of the amplitude M 0 is given by the ππ scattering phase δ. Indeed, employing the Sokhotski-Plemelj formula, where p.v. stands for the principal value, one can re-write Eq. (28)  where it was used that 1 + ie iδ sin δ = e iδ cos δ.
Since the expression in parentheses is real by construction, the entire phase of the amplitude M 0 is indeed given by the ππ scattering phase δ. Generalization of the formalism to the coupled-channel case is straightforward and amounts to replacing the scalar objects by matrices, as given in Eq. (29) in Sec. III A. σ π (s) Therefore, studies of the amplitude M L 0 (s) amount to building an analytical form of the discontinuity Disc C 0 (s) which reproduces the tabulated triangle loop function C 0 (s) in the kinematical regimes of interest through the dispersive integral (B2). Υ Υ ′ π π Z b π π Figure 6. The triangle diagram with the π-π FSI described by the scalar loop function C0 as introduced in Eq. (B1). The dotted line indicates the cut -see Eq. (B4).