Spin partners $W_{bJ}$ from the line shapes of the $Z_b(10610)$ and $Z_b(10650)$

In a recent paper Phys.Rev. D98, 074023 (2018), the most up-to-date experimental data for all measured production and decay channels of the bottomonium-like states $Z_b(10610)$ and $Z_b(10650)$ were analysed in a field-theoretical coupled-channel approach which respects analyticity and unitarity and incorporates both the pion exchange as well as a short-ranged potential nonperturbatively. All parameters of the interaction were fixed directly from data, and pole positions for both $Z_b$ states were determined. In this work we employ the same approach to predict in a parameter-free way the pole positions and the line shapes in the elastic and inelastic channels of the (still to be discovered) spin partners of the $Z_b$ states. They are conventionally referred to as $W_{bJ}$'s with the quantum numbers $J^{PC}=J^{++}$ ($J=0,1,2$). It is demonstrated that the results of our most advanced pionful fit, which gives the best $\chi^2/{\rm d.o.f.}$ for the data in the $Z_b$ channels, are consistent with all $W_{bJ}$ states being above-threshold resonances which manifest themselves as well pronounced hump structures in the line shapes. On the contrary, in the pionless approach, all $W_{bJ}$'s are virtual states which can be seen as enhanced threshold cusps in the inelastic line shapes. Since the two above scenarios provide different imprints on the observables, the role of the one-pion exchange in the $B^{(*)}\bar{B}^{(*)}$ systems can be inferred from the once available experimental data directly.


I. INTRODUCTION
Heavy-quark spin symmetry (HQSS) is an approximate symmetry of QCD. It states that in the limit of an infinite mass of the heavy quark its spin is conserved by the strong interactions. As one departs from this limit, corrections scale as Λ QCD =M Q , where Λ QCD ≃ 200 MeV denotes the intrinsic mass scale of QCD and M Q is the mass of the heavy quark. Stated differently, in the heavy-quark limit hadronic interactions do not depend on the heavy-quark spin orientation, and hadronic states can be classified by the quantum numbers of the light degrees of freedom (d.o.f.). This implies that an experimentally observed state with a given heavy quark spin should have spin partner states with different heavy quark spins but the same light quark cloud. This pattern of spin partners is well established amongst the states with the assumedQQ structure below the open-flavor threshold. For example, in the b-sector, for each radial excitation number n, one identifies η b ðnSÞ (with the heavy quark-antiquark spin S QQ ¼ 0) as a spin partner of the ΥðnSÞ (S QQ ¼ 1), and h b ðnPÞ (S QQ ¼ 0), as a spin partner of the χ bJ ðnPÞ (J ¼ 0; 1; 2) (S QQ ¼ 1). Meanwhile, the present experimental situation does not allow one to reliably identify such spin partners for the confirmed exotic states in the b-sector above the open-flavor threshold, namely, for the Z b ð10610Þ and Z b ð10650Þ states. However, the high luminosity and high statistics Belle-II experiment which has just started to operate [1] should be able to provide us with new information on these exciting states.
This work contains a set of model-independent predictions based on HQSS for both pole parameters and line shapes of the charged J PC ¼ J þþ (J ¼ 0, 1, 2) molecular states in the spectrum of bottomonium. The parameters of the formalism were fixed in an earlier analysis of the most recent experimental information available for the 1 þ− bottomoniumlike states Z AE b ð10610Þ and Z AE b ð10650Þ (for simplicity often referred to as Z b and Z 0 b , respectively). They were observed by the Belle Collaboration as peaks in the invariant mass distributions of the ϒðnSÞπ AE (n ¼ 1, 2, 3) and h b ðmPÞπ AE (m ¼ 1, 2) subsystems in the dipion transitions from the vector bottomonium ϒð10860Þ [2] and later confirmed in the elastic B ðÃÞBÃ channels [3][4][5]. Being isovectors, the Z AE b ð10610Þ and Z AE b ð10650Þ cannot be conventionalbb mesons as their minimal possible quark content is four-quark while their proximity to the BB Ã and B ÃBÃ thresholds, respectively, provides a strong support for their molecular interpretation-see, e.g., recent review articles [6,7]. In particular, the interference pattern in their inelastic decay channels ϒðnSÞπ and h b ðmPÞπ can be explained naturally in the framework of the molecular picture [8]. For a competing tetraquark interpretation claimed to be also compatible with the data see Refs. [9][10][11].
Since the mass of the b quark is very large compared with the typical QCD scale, M b ≫ Λ QCD , the constraints from the heavy-quark spin symmetry should be very accurate for bottomoniumlike systems, including the Z AE b ð10610Þ and Z AE b ð10650Þ. The wave functions of the Z b states in the molecular picture can be written as [8] where S P qq denotes the wave function of the light qq pair with the total spin S and parity P, and S P bb means the same for the bb pair. Based on this structure four additional isovector sibling states W bJ are predicted to exist [7,8,[12][13][14] with the quantum numbers J PC ¼ J þþ (J ¼ 0, 1, 2) and with the wave functions Thus, using that the low-energy interaction and the transition potentials between the elastic and inelastic channels in the Z b 's fixed from the existing experimental data can be uniquely translated into the W bJ sector using HQSS constraints, the theoretical description of the spin partner states W bJ appears to be rather straightforward. In Ref. [15] an effective field theory (EFT) approach to the Z b ð10610Þ and Z b ð10650Þ was developed and employed in the analysis of the experimental data on the line shapes of these states in the elastic and inelastic channels. The approach is formulated based on an effective Lagrangian consistent with both chiral and heavy-quark spin symmetry of QCD. The key features of the approach and the central findings of Ref. [15] can be summarized as follows (see Ref. [16] for a review of a similar chiral EFT approach in few-nucleon systems): (i) The EFT is constructed employing the so-called Weinberg counting [17,18], proposed originally to treat few-nucleon systems. The potential is constructed to a given order in Q=Λ h (here Q denotes the soft scales of the given problem and Λ h ≈ 1 GeV represents the hard scale of the chiral EFT) and then resummed nonperturbatively employing the Lippmann-Schwinger equation. Accordingly, at leading order (LO) the potential contains two momentumindependent, OðQ 0 Þ, contact interactions and the pion exchange. (ii) Simultaneously with the chiral EFT expansion, the potential is expanded around the spin symmetry limit. At leading order in chiral EFT this calls for the inclusion of the B Ã -B mass difference together with all interaction vertices constructed in line with HQSS. Since Λ QCD =M b ≈ 0.04 ≪ 1, subleading HQSS violating contributions, which would lead to additional terms in the potential, are expected to play a minor role. This expectation was confirmed in Ref. [15], where it was shown that the data in the 1 þ− channel were essentially consistent with HQSS constraints imposed on the potential. (iii) The binding momenta, the pion mass and the momentum scale generated by the splitting between the BB Ã and B ÃBÃ thresholds [δ ¼ m Ã − m ≈ 45 MeV with m (m Ã ) denoting the B (B Ã ) meson mass], are treated as soft scales of the system, generically called Q above. (iv) In order to remove the strong regulator dependence caused by the high-momentum contributions from the S-wave-to-D-wave B ðÃÞBðÃÞ → B ðÃÞBðÃÞ transitions (in what follows, simply S-D transitions) generated by the one-pion exchange (OPE) a promotion of the OðQ 2 Þ S-wave-to-D-wave counterterm to leading order is required. Then, the fit to the data enforces that a large portion of the S-D contribution from the OPE gets balanced by the S-D contact interaction. However, the residual effect from the OPE on the line shapes is visible and results in a quantitative improvement of the fits. To check the convergence of the scheme the effect from the other contact interactions at the order OðQ 2 Þ, namely, from two S-wave-to-S-wave terms, was also studied. However, their effect on the line shapes was shown to be numerically small in line with the assumed power counting. (v) The effect of the inelastic channels ϒðnSÞπ (n ¼ 1, 2, 3) and h b ðmPÞπ (m ¼ 1, 2) is included by allowing them to couple to the S-wave B ðÃÞ -meson pairs. Following Refs. [19,20], transitions between inelastic channels are omitted in the potential. As a result, the effective elastic potentials acquire imaginary parts driven by unitarity while the contributions to the real parts of the elastic potentials can be absorbed into the redefinition of the momentumindependent OðQ 0 Þ contact interactions. The treatment of the inelastic channels used in this work is analogous to the construction of the annihilation potential in nucleon-antinucleon scattering-see, e.g., Ref. [21]. (vi) Extension of the approach to the SU(3) sector requires that all the other members of the lightest pseudoscalar Goldstone-boson octet be treated also as explicit d.o.f. already at leading order. That is why the one-η exchange (OEE) is also included as a part of the B ðÃÞBðÃÞ → B ðÃÞBðÃÞ effective potential. In the SU(2) sector, however, the effect from the explicit treatment of the η-meson appears to be negligible. (vii) All low-energy constants (the two elastic couplings, the effective couplings to the inelastic channels, as well as the S-S and S-D contact interactions) are fixed from a combined fit to the experimental line shapes in the decays ϒð10860Þ → BB Ã π, B ÃBÃ π, h b ð1PÞππ, and h b ð2PÞππ which proceed via the excitation of the Z b ð10610Þ and Z b ð10650Þ exotic states as well as from the total rates for the decays ϒð10860Þ → ϒðnSÞππðn ¼ 1;2;3Þ. The line shapes in the ϒð10860Þ → ϒðnSÞππ channels could not be included in the analysis so far since they require a proper treatment of the two-pion final-state interaction. It is instructive to comment on the difference between the EFT used in this work and the heavy-quark effective theory (HQET) (see, e.g., the textbook [22]). In the HQET framework, single heavy-quark (heavy-light) systems are considered and the d.o.f. related to the heavy quark can be integrated out explicitly which makes the velocity power counting of operators manifest. Meanwhile, the power counting in the systems containing two heavy hadrons, like in the B ðÃÞBðÃÞ systems considered in this work, is quite nontrivial and substantially different from that in the processes involving a single heavy particle. The reason for that is the presence of pinch singularities which show up in loop contributions containing a B ðÃÞBðÃÞ cut, as soon as the heavy particles are treated as static sources [23]. The problem of the pinch singularity is identical to the one discussed extensively in the context of chiral EFTs for the two-nucleon system [18]. To cure this problem, the term ∇ 2 =ð2MÞ has to be included in leading-order calculations of the observables (together with ∂ 0 ), as shown by Weinberg in Ref. [18]. A consequence of this different behavior of heavy-heavy molecular systems is that one cannot relate different heavy-quark sectors within the same (heavy-flavor) EFT [24]. As a result, the findings reported in this work for the spin partners in the b-quark sector cannot be translated to the c-quark sector in a controllable way.
It is also important to emphasize that, in our EFT framework aiming at a model-independent understanding of the near-threshold states in a limited energy domain (of the order of δ around the BB Ã threshold), we cannot probe the details of various short-range interactions. Therefore, we do not include compact quark components of the wave functions explicitly but take them into account (as well as other possible short-range contributions) effectively via a series of contact interactions with different numbers of derivatives. Thus, in the EFT approach, a bare pole ∝ 1=ðE − E 0 Þ is represented by a series in powers of ðE=E 0 Þ n ðn ¼ 0; 1; …Þ, which is then absorbed by the tower of contact interactions. Therefore, if a bare pole were necessary within the regime of applicability of the EFT, the effective potential would strongly depend on the energy, which would imply that the EFT and in particular the derivative expansion for the contact terms would not converge. However, this is not the case for the bottomoniumlike states Z b ð10610Þ and Z b ð10650Þ. Indeed, it follows from the data analysis performed in Ref. [15] (and from the discussion below) that the existing experimental data are well described without the necessity of including compact quark components explicitly. This could be anticipated since the inclusion of the S-wave-to-S-wave contact terms with two-derivatives (which enter our EFT at next-to-leading (NLO) order) indeed results in a perturbative effect, in line with our power counting, and, therefore, indicates convergence of the EFT series.
In this work, we provide parameter-free predictions for the HQSS partner states of the Z b and Z 0 b molecules using the framework summarized above. In particular, we exploit the fact that all the parameters of the elastic and inelastic potentials extracted from the experimental line shapes in the J PC ¼ 1 þ− channel are the same in the partner channels up to spin symmetry violating corrections in the contact interactions that are expected to be small. Thus we are able to predict, for the first time, the line shapes and to extract the pole positions for the spin partner states W bJ with the quantum numbers J þþ (J ¼ 0, 1, 2). In particular, we discuss the impact of the one-pion exchange on the observables.
It needs to be mentioned that the same formalism for the spin partner states of the Z ð0Þ b 's was employed in Ref. [14]. However, compared with that paper, this work marks progress in four important aspects: (i) All parameters of the interaction are now fixed directly from a fit to the measured line shapes contrary to the earlier study where the masses of the Z b states obtained in different analyses were used as input. (ii) Relevant inelastic channels are included which makes it possible to predict the line shapes in the inelastic channels, in addition to the elastic ones. (iii) We investigate how the renormalization program works in the coupled-channel case. In particular, we are now in a position to study the effect of the S-D counterterm on the line shapes and the pole locations of the spin partner states. The need for this counterterm is one of the conclusions of Ref. [15]. (iv) The uncertainty of the EFT predictions for the pole positions, which comes from various sources, is estimated and discussed. The paper is organized as follows. Section I contains a brief introduction to the formalism employed. In Sec. II, the effective potentials in the Z b 's and W bJ 's channels are discussed in detail. In Sec. III, the coupled-channel equations are provided and the expressions for the differential production rates in all elastic and inelastic channels are constructed. The resulting formulas are then employed in Sec. IV to analyze the existing experimental data in the Z b 's channel (in line with the results of Ref. [15]) and to predict the line shapes in the W bJ 's channels with the quantum numbers J þþ (J ¼ 0, 1, 2). In addition, in Sec. V, the pole parameters (locations and residues) for the W bJ states are provided and an analysis of uncertainties is presented. We summarize in Sec. VI. Appendix A contains the details of the NLO Lagrangian OðQ 2 Þ used to build the suitable EFT, while Appendix B provides the details of the partial wave projection operators applied to the effective potential.

A. Some generalities and definitions
The partial-wave-projected effective potential in the elastic channels derived in the EFT framework and used in our calculations reads where V CT eff , V π and V η stand for the effective contact interaction potential (composed of the elastic, V CT NLO , and inelastic, δV, contributions, as discussed below), OPE, and OEE, respectively, and the indices α and β, which depend on the particle channel and quantum numbers (J PC ), are defined as Here the individual partial waves are labeled as 2Sþ1 L J with S, L, and J denoting the total spin, the angular momentum, and the total momentum of the two-meson system, respectively. Finally, the sign in the parentheses corresponds to the BB Ã states with a given C-parity with a universal definition of the C-parity transformation employed,Ĉ for any meson M.

B. Contact interactions
The OðQ 0 Þ short-ranged elastic (open-bottom) interaction between the states with given quantum numbers composed of the B ðÃÞBðÃÞ pairs is parametrized in terms of two contact terms C 10 and C 11 ; see Refs. [25,26] and Lagrangian (A1) quoted in Appendix A. This appendix also contains the momentum-dependent order OðQ 2 Þ contact interactions shown by Eq. (A2) and originally derived in Ref. [15]. Formally, the order OðQ 2 Þ chiral Lagrangian contains also the contact terms which scale with the pion mass squared (m 2 π ). However, as long as we work at a fixed light-quark mass those can be absorbed into the leadingorder counterterms. Thus, to order OðQ 2 Þ the contact potentials in the elastic channels for various quantum numbers relevant for this study read where each potential is given for the basis states defined by Eq. (9) and p (p 0 ) stands for the relative momentum of the initial (final) heavy-meson pair. The potential given in Eq. (12) was derived and used already in Ref. [15].

C. Inelastic channels
Since the interaction between the pion and heavy quarkonia is suppressed (see the discussion in Refs. [19,20]), the direct transitions between the inelastic (hidden-bottom) channels can be safely neglected in the potential. Based on this assumption, in the approach employed in Refs. [15,19,20], the effect of the inelastic channels on observables is included through the transitions between the inelastic and elastic potentials only, while unitarity is preserved. Furthermore, it is argued in Ref. [15] that all inelastic channels only couple to the S-wave elastic ones as their couplings to the D-wave elastic channels are suppressed by the factor p 2 typ =m 2 ≪ 1. Transitions between the S-wave elastic and inelastic channels are described by the Lagrangian [13] L inel HH ¼ Here the spin multiplets of the heavy-light mesons read where σ i are the Pauli matrices, P a ðP a Þ and V i a ðV i a Þ are the pseudoscalar B (B) and vector B Ã (B Ã ) mesons, respectively, with a and b for the isospin indices. The multiplets of the heavy ðbbÞ mesons are built as and with f π ¼ 92.4 MeV being the pion decay constant [27]. Since the quarkonium states (ϒ,η b ) and (h b , χ b ) form spin multiplets [see Eqs. (18) and (19)], the coupling constants fixed in the analysis of the line shapes in the 1 þ− channels [15] can be used to account for the inelastic transitions in the spin partner channels. As long as the direct inelastic transitions are neglected, the effect of the inelastic channels on the elastic ones can be included via an additional contribution, δV, to the effective contact elastic-to-elastic transition potential [19], where Here P αβ stands for the real part of δV αβ , m h i denotes the mass of the heavybb meson in the ith inelastic channel and M is the total energy of the system. Further, v iα ðp i ; pÞ is the vertex function for the transitions between various heavy-meson states (9) (labeled by greek letters α, β and so on) and inelastic channels (labeled by latin letters i, j and so on). The arguments p i and p denote the on-shell momenta of the inelastic and elastic channels involved, respectively, measured in the rest frame of the system. One finds for i ¼ ϒðnSÞ, where E π ðp i Þ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi m 2 π þ p 2 i p denotes the pion energy for a given inelastic momentum, with λðx; y; zÞ being the standard triangle function.
In contrast to Ref. [15], we now keep the energy dependence in the vertices v ϒ iα explicitly, as it comes out from the Lagrangian (16). However, since the variation of the inelastic momenta with the energy is very minor near the elastic thresholds, this correction does not affect the quality of the fits and merely results in rescaling of the inelastic coupling constants compared with those used in Ref. [15]. For i ¼ h b ðmPÞ, χ bJ ðmPÞ, with m ¼ 1, 2 and J ¼ 0, 1, 2, the expression for the vertex reads The coefficients ξ ϒ iα and ξ χ iα , provided explicitly in Table I, are fixed by the HQSS and are straightforwardly calculated from the traces appearing in Eq. (16). They do not depend on n and m although the individual coupling constants, in principle, do. We also note that the relative signs between various couplings are only relevant in the particle coupled channels (the channels with J PC ¼ 1 þ− and 0 þþ in Table I while for the channels J PC ¼ 1 þþ and 2 þþ only the absolute values of ξ's enter).
The real parts induced by the inelastic channels, P αβ , are divergent and need to be regularized. The scheme employed in Ref. [15] assumes that the whole real part of the inelastic contribution in the J PC ¼ 1 þ− channel is absorbed into a redefinition of the low-energy constants (LECs) C d and C f ; see Eq. (12). This is justified as the momentum dependence of P αβ coming from remote inelastic channels is very weak and, therefore, can be neglected. To proceed we need to ensure that, in the heavyquark limit, the same approach works for the complete spin multiplet. Then we have TABLE I. HQSS-constrained coefficients ξ ϒ iα and ξ χ iα in front of the coupling constants in the multiplets ϒ and χ calculated explicitly from the traces in Eq. (16)-see the vertices given in Eqs. (24) and (26). The symbol Á Á Á indicates that the corresponding channels do not couple with each other.

Multiplet
Channel These principal value integrals are factored out of the brackets in Eq. (28) since the masses of the members of the spin multiplets coincide in the heavy-quark limit. One finds by a direct evaluation using the coefficients from Table I (no summation in α or β is implied here) for α ¼ BBð 1 S 0 Þ and β ¼ B ÃBÃ ð 1 S 0 Þ; and finally where α ¼ B ÃBÃ ð 5 S 2 Þ. For the 1 þþ and 2 þþ channels, the individual contributions from the intermediate states πχ bJ with J ¼ 0, 1, 2, in order, are quoted explicitly on the righthand side (rhs) of Eqs. (32) and (33). According to Eqs. (22), (23), (28) and (30) the real parts of the loops in the 1 þ− channel can be absorbed into the bare LECs C d and C f entering the short-range interaction for the It is also straightforward to see using Eqs. (28) and (31)-(33) that these redefinitions also hold for the interactions in the spin partner channels J þþ (J ¼ 0, 1, 2). This procedure is correct up to the neglected terms that violate spin symmetry and contain the energy dependence of the integrals in Eq. (29). Thus, in line with Ref. [15], in what follows only the imaginary parts of the inelastic loops are retained in the effective contact interaction potential (22) for all spin partner states.

D. Pion exchange
The pion exchange in the B ðÃÞ B Ã system is described by the Lagrangian [28,29] where A ¼ i ffiffi 2 p ðϵ × ϵ 0Ã Þ, with ϵ and ϵ 0Ã the polarization vectors of the initial and final B Ã mesons, and q is the pion momentum. These vertices agree with those used in Ref. [30].
In order to determine the dimensionless coupling constant g b we rely on the heavy-quark flavor symmetry and set where the numerical value of the g c is extracted from the most recent measurement of the D Ãþ → D 0 π þ decay width, Here m D Ãc and m D 0 denote the masses of the D Ãþ and D 0 mesons, respectively, and the final-state momentum q ¼ 39 MeV [27]. The value of g b quoted in Eq. (37) agrees within 10% with the result of a recent lattice QCD determination of the B Ã Bπ coupling constant [31]. The isospin factor for the OPE potential is τ 1 · τ c 2 ¼ −τ 1 · τ 2 ¼ 3 − 2IðI þ 1Þ that gives (−1) for the isotriplet states considered in this work. Here τ c ¼ τ 2 τ T τ 2 ¼ −τ is the charge-conjugated Pauli matrix used for the antifundamental representation of the isospin group. Finally, the overall sign of the OPE potential depends on the C-parity of the channel-see, e.g., Ref. [32] for a detailed discussion. Using the definition of the C-parity given in Eq. (10) one finds where V π ≡ hBB Ã jV π jBB Ã i ¼ hBB Ã jV π jBB Ã i: ð40Þ We consider a coupled-channel system for the BB, BB Ã =BB Ã and B ÃBÃ channels. Using the labels one can write for the OPE potentials where the contributions from both time orderings as obtained in time-ordered perturbation theory (TOPT) are taken into account (see Figs. 1 and 2). Further, q ¼ p − p 0 , ϵ 1 and ϵ 2 (ϵ 0 1 Ã and ϵ 0 2 Ã ) stand for the polarization vectors of the initial (final) B Ã mesons, The denominators D B ðÃÞ B ðÃÞ π ðp; p 0 Þ correspond to the B ðÃÞ B ðÃÞ π propagators written in TOPT for the nonrelativistic B and B Ã mesons, The time-reversed transition potentials V π 31 ðp; p 0 Þ and V π 32 ðp; p 0 Þ are trivially obtained from Eqs. (43) and (45) by interchanging the particle labels as 1 ↔ 1 0 , 2 ↔ 2 0 (see also Fig. 2).
Since our analysis covers the energy region between the BB and B ÃBÃ thresholds split by approximately 2δ ≈ 90 MeV, which numerically appears to be of the order of the pion mass, a relativistic expression for the pion energy is used. It should be noted that, in the region of interest around the B ðÃÞBðÃÞ thresholds, M < 2m þ m π , so that, unlike charmonium systems, one never hits the threebody cut in the OPE potentials defined above.
Since the pion is emitted by B ðÃÞ mesons in the P-wavesee the Lagrangian (34) and the vertices (35) and (36)-the OPE potential mixes S and D waves. The partial wave projection of the OPE potentials (43)- (46) can be done using the formalism of Refs. [14,33,34] that gives where L ¼ S, D and G, n ¼ p=p (n 0 ¼ p 0 =p 0 ), and a complete set of relevant properly normalized projection operators PðJLS; nÞ is given in Appendix B. Note also that, as was stressed in Ref. [35], the OPE (as well as OEE) potential in the system of two heavy-light pseudoscalar and vector mesons is well defined in the sense of an effective field theory only in connection with the contact operator provided by the short-range potential V CT eff [Eqs. (12)-(15)].

A. Production vertex
The twin states Z b ð10610Þ and Z b ð10650Þ are produced in the one-pion decays of the ϒð10860Þ resonance, Because of a different C-parity, the spin partner states W bJ should be produced in radiative decays of the ϒð10860Þ, ϒð10860Þ → γW bJ → final state: Production and decay channels for the Z b 's and W bJ 's taken into account in our approach are summarized in Fig. 3. In line with the discussion in Sec. II C, since the couplings of the ϒð10860Þγ source term with the D-wave elastic channels are suppressed, we retain only the couplings of the ϒð10860Þγ source term with the elastic channels in the S wave.
A set of diagrams which contribute to the process ϒð10860Þ → γB ðÃÞ B ðÃÞ and provide a gauge invariant amplitude is shown in Fig. 4, where diagrams (a), (b1) and (b2) contribute to the production operators at tree level while diagrams (c)-(e) represent contributions from the loops assuming that the intermediate particles are B and B Ã mesons only. In this work, we do not aim at predicting the absolute rate of the decays ϒð10860Þ → γW bJ , which might involve some more sophisticated mechanisms (e.g., as advocated in Ref. [36]), but rather focus on the energy dependence of the line shapes with an arbitrary overall normalization. To this end, we assume that for this process (similarly to the Z b 's production) the energy dependence from the production operator is rather smooth in the vicinity of thresholds and can be merely neglected as compared to rapidly varying B-meson amplitudes in the final state. To advocate this approximation, below we discuss the diagrams shown in Fig. 4 in more detail.
To estimate the strength and the structure of the source term we start from the effective Lagrangian connecting the ϒð10860Þ bottomonium with the heavy-meson fields at leading order in the HQSS expansion [13], FIG. 3. Summary of the production and decay channels for the Z b 's and their spin partners W bJ 's considered in this work. The states and thresholds are arranged from bottom to top in accordance with the increasing energy.
Gauging this Lagrangian leads to a set of contact ϒð10860Þ → γB ðÃÞBðÃÞ vertices which contribute to diagrams (a) and (d), where, for a given J, index α runs over the relevant S-wave states, that is, The spin operators, normalized according to P λ jŜ ðJ þþ Þ α Here ϵ γÃ , ϵ ϒ and ϵ 1ð2Þ denote the polarization vectors of the photon, ϒ and B Ã mesons, respectively, and the explicit forms of the projectors P on relevant heavy-meson states are given in Appendix B. Further, the partial-wave-projected vertices v where e is the magnitude of the electron charge and the ratios of the coupling constants, λ ðJ þþ Þ α , related by HQSS are quoted in Table II. It is shown in Ref. [37] that experimental data might call for a significant amount of spin symmetry violation in the transition ϒð10860Þ → B ðÃÞBðÃÞ (there is a tension of 2σ between the spin symmetric ratio of decay widths and the experimental data). Since the same couplings also contribute to the transitions ϒð10860Þ → γB ðÃÞBðÃÞ , HQSS violation is a potential additional source of uncertainty for our results-we come back to this issue in the discussion below.
In addition to the contact diagram with the photon emission from the ϒð10860Þ → γB ðÃÞBðÃÞ vertex, diagrams with the photon emission from the B ðÃÞ -meson lines should be considered with the B ðÃÞ → B ðÃÞ γ vertices being of an electric or magnetic type. While the amplitudes with the magnetic photon emission vertices are gauge invariant by themselves, the additional amplitudes with the electric photon emission vertices are important to compensate for the gauge dependence of the contact W bJ → γB ðÃÞBðÃÞ diagram and thus to provide an overall gauge invariance of the full amplitude. To estimate the electric contributions, we notice that in the nonrelativistic heavy-meson formalism used here all momenta involved are 3-momenta and the photon momentum k fulfills the relation ϵ γ · k ¼ 0. Then, one readily arrives at the following estimates for the treelevel diagrams (b1) and (b2) with the electric (thence superscript e) photon emission from the external B ðÃÞ -meson lines relative to the contact W bJ → γB ðÃÞBðÃÞ diagram, where one power of the relative B ðÃÞ -meson momentum p α comes from the ϒð10860Þ → B ðÃÞBðÃÞ vertex extracted from Eq. (51); the term 1=ðmωÞ is the B ðÃÞ -meson propagator with ω for the photon energy (here we do not distinguish between the B and B Ã mass); and the second factor p α comes from the  , responsible for the production of the W bJ states in the radiative decays ϒð10860Þ → γW bJ .
20=3 p electric photon emission vertex from the B ðÃÞ meson, which is derived by gauging the kinetic term in the Lagrangian (A1). Further, to arrive at the very last relation in Eq. (58) we used that the energy of the B-meson pair relative to the threshold in the channel α, E α ≈ p 2 α =m, does not exceed several dozen MeV while the photon energy ω is an order of magnitude larger, ω ≈ M ϒð10860Þ − 2m ≈ 200-300 MeV.
One is, therefore, led to conclude that the diagrams (b1) and (b2) (with the electric photon emission from the B ðÃÞmeson lines) while being important to guarantee gauge invariance of the amplitude, in practical calculations provide only small corrections. Thus, the tree-level amplitude behaves basically as a constant in the energy region of relevance.
The loop contributions from the diagrams (c)-(e) were already studied in the literature in the context of scalar mesons made of light quarks-see, e.g., Refs. [38,39]. In particular, it is shown that for pseudoscalar mesons such loops form a gauge invariant subset of diagrams which yields a finite contribution to the amplitude. The arguments of Refs. [38,39] can be generalized to find that, in the HQSS limit, these conclusions hold also for all members of the heavy-meson spin multiplet and for all quantum numbers J þþ . Further, using the explicit results of Refs. [38,39] for the loops with a pointlike interaction between the mesons in the final state one concludes that, to a good approximation, also for diagrams (c)-(e) the production operator can be treated as a constant.
To illustrate the argument, consider the resulting contribution from the diagrams shown in Fig. 4 for the uncoupled case, where the vertex and the spin structure in front of the parenthesis are from Eqs. (55)-(57), A α ðp α Þ and ip α denote the real and imaginary parts of the pertinent loop, and f α on ðp α Þ is the on-shell B-meson amplitude in the final state. Unitarity forces f α on to have the form where B α ðp α Þ denotes the real part of the inverse scattering amplitude which is unconstrained by unitarity and is a real meromorphic function of p 2 α near the origin p α ¼ 0. To leading order in a momentum expansion, B α ðp α Þ is given by the inverse scattering length. Then one finds Thus, unitarity forces the production amplitude to be proportional to the scattering amplitude in the final state (a coupled-channel version of this relation is provided in Ref. [40]). In the heavy-quark limit the functions A and B do not depend on the channel. Moreover, since the momentum dependence of the functions Aðp α Þ and B α ðp α Þ is controlled by the left-hand cuts of the production operator and the scattering amplitude, respectively, we expect that near thresholds both are well approximated by constants, which are also independent of the channel in the heavy-quark limit. Based on this one can predict the ratios of the partial widths for different decay channels of the W bJ 's, up to spin symmetry violating corrections. It is proposed in Ref. [36] that the most prominent production mechanism for the Z b states in the ϒð10860Þ and ϒð11020Þ decays involves B 0 1B or B 0B intermediate states, with B 0 and B 0 1 being the broad members of the quadruplet of the positive P-parity B mesons. If this proposal is correct, the decay mechanism through the B ðÃÞBðÃÞ pairs considered above will give only a small contribution. However, it should be stressed that the mechanism proposed in Ref. [36] should not change the line shapes but only the total rate of the production cross sections, which is not a subject of the current study.

B. Coupled-channel system
The set of the allowed quantum numbers for the B ðÃÞBðÃÞ system is encoded in the basis vectors quoted in Eq. (9). Inclusion of the OPE interaction enables transitions to the D and even G waves [33].
For a given set J PC the system of the partial-wavedecomposed coupled-channel Lippmann-Schwinger-type equations reads where α, β, and γ label the basis vectors defined in Eq. (9); the effective potential is defined by Eq. (8); and the scattering amplitude T αβ is related with the invariant amplitude M αβ as with m 1;α and m 2;α (m 1;β and m 2;β ) being the masses of the B ðÃÞ mesons in the channel α (β). The two-body propagator for the given set J PC takes the form where the reduced mass is It is convenient to define the energy E i relative to a particular threshold, namely, Finally, to render the loop integrals well defined we introduce a sharp ultraviolet cutoff Λ which needs to be larger than all typical three-momenta related to the coupledchannel dynamics. For the results presented below we choose Λ ¼ 1 GeV but we also address the problem of the renormalizability of the resulting EFT and estimate and discuss the theoretical uncertainty from the cutoff variation.

C. Production rates
Since we are not interested in the absolute scale but only in the energy dependence of the line shapes, the production amplitude of the βth elastic channel from a pointlike source for some given quantum numbers J PC can be defined as where the relativistic normalization factor is and the nonvanishing partial-wave-projected production vertices v α are from Eq. (57).
Since the direct interactions between the inelastic channels are neglected in the formalism applied here, the ith inelastic channel in the final state can only be reached via a transition through the intermediate elastic channels. In particular, for a given set J PC , the inelastic amplitude M i is obtained by convolving the relevant elastic amplitude M β ðM; pÞ from Eq. (67) with the corresponding elastic-to-inelastic transition vertex from Eqs. (24) and (26), that is, where Finally, the differential widths in the elastic and inelastic channels read respectively, where k is the three-momentum of the photon in the rest frame of the ϒð10860Þ and p β is the threemomentum in the βth elastic channel in the rest frame of the B ÃBðÃÞ system, namely, In what follows, we consider three fitting strategies introduced in Ref. [15]: (1) Contact fit: Purely S-wave momentum-independent contact interactions (analogous to fit A in Ref. [15]). (2) Pionful fit 1: Complete leading-order potential that involves S-wave contact terms plus the OPE, plus the OðQ 2 Þ S-wave-to-D-wave counterterm promoted to LO (analogous to fit E in Ref. [15]). (3) Pionful fit 2: Pionful fit 1 supplemented by the OðQ 2 Þ S-wave-to-S-wave contact terms at NLO and the η-meson exchange (analogous to fit G in Ref. [15]). The line shapes in the 1 þ− channel, where the Z b 's states reside, corresponding to the best fits for the three schemes quoted above, are compared with the experimental data in Fig. 5. The parameters extracted from these fits are collected in Table III. One can see that the quality of the line shape description by the pionful fits is better than that by the contact fit that is reflected in the change of the χ 2 =d:o:f: from 1.29 for the contact fit to 0.95 for the pionful fit 1 and 0.83 for the pionful fit 2.

B. Renormalizability of the heavy-hadron EFT with pions
The use of the standard nonrelativistic approach to heavy mesons, as employed in Ref. [15] and also used here, leads to coupled-channel integral equations for the scattering amplitudes which, at leading order in the EFT expansion, are linearly divergent. Therefore, when the potential truncated at a given order is iterated within the integral equations an infinite number of ultraviolet (UV) divergent higher-order contributions is generated. The problem is well known in the context of nuclear chiral EFT-see, e.g., Refs. [41,42] and references therein. The standard way to cure this problem in practical calculations is to employ a finite UV cutoff of the order of a natural hard scale in the problem, so that the unwanted higher-order contributions turn out to be suppressed [43]; see also a recent discussion in Ref. [44]. For an alternative approach with relativized integral equations of the Kadyshevsky type in the context of a nucleon-nucleon (NN) and heavy-meson EFTs see Refs. [45] and [46], respectively.
The logic explained above is the basis for the renormalization program used in Ref. [15] and, hence, is also employed here. It should be stressed that the formulation of an EFT for the Z b 's and their spin partners is much more challenging than that for the NN problem because of the larger soft scales involved here. Indeed, since the B mesons are, roughly, by a factor of 5 heavier than nucleons, an EFT for the Z b ð10610Þ and Z b ð10650Þ states, separated by δ ¼ 45 MeV, unavoidably involves the momenta of the order of ffiffiffiffiffiffi mδ p ≈ 500 MeV treated as soft. Moreover, in the course of practical fits of the experimental line shapes the momenta as large as 800 MeV are included from the high-energy tail of the experimental distributions. Clearly, the influence of such high momenta on the dynamics close to the relevant thresholds is minor whereas the renormalization of the theory (and the possible residual cutoff dependence) is severely affected by this high-momentum range. In Ref. [15] it was found that the cutoff dependence generated via the iterations of the OPE in the S waves can be almost completely absorbed into the momentumindependent contact terms C d and C f at LO. On the other hand, the cutoff dependence of the line shapes in the 1 þ− channel from iterations of the S-wave-to-D-wave OPE turns out to be sizable which calls for the promotion to leading order of the contact term D SD providing the S-wave-to-D-wave transitions between heavy mesons, which naïvely would appear only at NLO. In Fig. 6, we illustrate the cutoff dependence for the elastic line shapes corresponding to the quantum numbers 1 þ− , 1 þþ and 0 þþ for the cutoff variation from 0.8 to 1.3 GeV. We start the discussion with the results for the quantum numbers 1 þ− , where the fits to the experimental data were performed (see the left plots in Fig. 6). As a general trend, the results demonstrate a mild cutoff dependence and a saturation for larger cutoffs. Meanwhile, as expected, the result for the smallest cutoff Λ ¼ 0.8 GeV for the pionful fit 1 deviates from the other curves (cf. the blue dashed curve for Λ ¼ 0.8 GeV with the red solid and dotted curves corresponding to Λ ¼ 1.0 GeV and Λ ¼ 1.2 GeV, respectively). Indeed, in order to maintain approximate Λ-independence in the pionful calculations for smaller cutoffs, we found empirically that the magnitude of the contact term D SD must be increased such that it generates an increasing S-wave-to-Swave higher-order contribution through iterations. The latter induces a strong Λ-dependence unless an additional OðQ 2 Þ S-S contact term is included in the potential. As a consequence, the results for the pionful fit 1 still show some cutoff dependence for the observables in the 1 þ− channel while the cutoff dependence for the pionful fit 2, where the order OðQ 2 Þ S-S contact term D d is included, is diminished significantly. Exactly the same pattern, though somewhat enhanced, can also be seen in Fig. 6 for the spin partners. It is obvious that, for the smallest cutoff, the results for the pionful fit 1 (blue dashed curve) possess an unwanted sizable S-wave-to-S-wave higher-order contribution which, for the pionful fit 2, is largely absorbed by the S-S contact term D d . Still, the results for the pionful fit 1 for the cutoffs from 1 GeVonward quickly saturate with the cutoff increase and may be regarded as reasonable predictions at leading order. In what follows, we will discuss the line shapes and extract the poles of the amplitude for both pionful fits 1 and 2. Still, we regard the predictions obtained for the pionful fit 2 as our main results since in this case the cutoff-related artifacts induced by the iteration of the truncated potential are significantly reduced. It should be clear that the results for the pionful fit 2 correspond to an incomplete NLO calculation and that, in addition, there are long-range contributions from the two-pion exchange (TPE) not included in the present study. It remains to be seen whether or not their inclusion affects the predictions for the spin partner states. However, given that, numerically, the longrange part of the OPE plays the role of a correction as compared with the short-range mechanisms, the effect from the long-range TPE is expected to be small.

C. Line shapes in the spin partner channels
In Figs. 7-9 the line shapes in the spin partner channels with J ¼ 0, 1, 2 are shown for the two pionful fit schemes introduced above. For each scheme the line shapes are calculated employing the best-fit parameters extracted from the analysis of the data in the 1 þ− channel for the cutoff 1 GeV. Specifically, in each plot we present the relevant elastic B ðÃÞBðÃÞ and inelastic differential widths (in arbitrary units) defined by Eq. (70). The relative normalization of the curves for the pionful fits 1 and 2 is chosen such that the two curves have the same magnitude at the resonance peak. In the case of the 0 þþ channel, where two states are present, the curves are normalized to have the same strength in one of the peaks. While the overall scale of the line shapes is not a subject of the current investigation, as discussed in Sec. III A, the branching fractions defined relative to the total width for each J are predicted here-see Table IV and discussions below. For the inelastic channels, as examples, we show the distributions in the χ b1 ð1PÞπ and η b0 ð1SÞπ final states. The energy behavior of the line shapes in the other channels not shown here [χ bJ ðmPÞπ with m ¼ 1, 2 and η b0 ð2SÞπ] is completely analogous to that for the χ b1 ð1PÞπ and η b0 ð1SÞπ channels, respectively, while their relative scales can be read off from Table IV. The differential rates for the pionful fits exhibit either a wellpronounced hump above the relevant threshold, as seen in Figs. 8 and 9 for the 1 þþ and 2 þþ channels, or a sizable near-threshold distortion, as in the 0 þþ case-see Fig. 7. The fitted values of the low-energy constants and couplings to the J PC ¼ 1 þ− data for the contact and pionful fits as defined at the beginning of Sec. IV. The OðQ 0 Þ contact terms C d and C f are given in units of GeV −2 , and the OðQ 2 Þ contact terms D D and D SD are in units of GeV −4 . The couplings g ϒðnSÞ (n ¼ 1, 2, 3) and g h b ðmPÞ (m ¼ 1, 2) are given in units of GeV −3=2 . Only the absolute values of the coupling constants are presented since physical quantities are not sensitive to their signs. Uncertainties correspond to a 1σ deviation in the parameters. The quality of each fit is assessed through the reduced χ 2 =d:o:f: quoted in the last column.
Fit This picture is typical for a resonance which is supported by the position of the poles of the amplitude in the energy complex plane-see the discussion in Sec. V below. In Fig. 10, for illustrative purposes, we compare the line shapes in the 2 þþ channel for the contact and pionful fit 2 schemes. For all J's, the inelastic line shapes corresponding to the contact fit reveal only a cusplike structure at the relevant elastic threshold enhanced by the presence of a near-threshold pole-the behavior typical for a virtual state scenario. On the contrary, when the pions supplemented by the Oðp 2 Þ contact terms are included, the poles move to the complex plane, as will be discussed in the next section, resulting in abovethreshold resonance-type structures in the line shapes.
The partial widths Γ in all considered elastic and inelastic channels can be obtained as integrals over the entire relevant energy interval. The ratios of the individual partial widths to the sum of all contributions for a given J are shown in Table IV. Such ratios do not depend on the overall scale and, therefore, can be regarded as a parameter-free prediction of our approach. In particular, for the elastic widths one finds the relations Although the potential spin symmetry violation in the elastic source terms discussed in Ref. [37] can somewhat distort these results, the general pattern should persist. Let us summarize the findings we arrived at.
(i) As expected in the molecular scenario, for each W bJ state, the decay rate to the corresponding elastic channel with the nearest threshold is the largest while the inelastic channels are strongly suppressed compared with it. The decay rates to remote elastic channels are also suppressed. For example, the contribution of the W 0 b0 state to the BB rate is quite marginal, as can be seen in Fig. 7. This is a direct consequence of the properties demonstrated by the data in the 1 þ− channel (see Fig. 5): although the FIG. 6. Propagation of the cutoff dependence of the theoretical fits from the 1 þ− channel used as input to the spin partner channels J þþ (J ¼ 0, 1, 2) which are parameter-free predictions. The upper panel shows the elastic line shapes for the pionful fit 1 (fit E in Ref. [15]) and the lower panel corresponds to the pionful fit 2 (fit G in Ref. [15]). In the upper panel, blue dashed curves denote Λ ¼ 0.8 GeV; red solid, Λ ¼ 1 GeV; red dotted, Λ ¼ 1. coupled-channel dynamics allows for such transitions, the data do not favor them. (ii) The largest rates correspond to the ϒð10860Þ decays to the γBB Ã and γB ÃBÃ channels via the W b1 and W b2 partners, respectively-see Eq. (73). (iii) The ratios predicted in Eq. (73) from the measured line shapes of the Z b states are consistent with the estimates presented in a recent study [47]. As to the absolute scale of the rates ϒð10860Þ→ γW bJ →γB ðÃÞBðÃÞ , a two-order-of-magnitude suppression is trivially expected as compared with the rates ϒð10860Þ → πZ ð0Þ b → πB ðÃÞBÃ because of the standard fine structure penalty for electromagnetic processes, which also agrees with the estimates made in Ref. [47]. Meanwhile, this suppression is expected to be overcome by the Belle-II experiment due to its large luminosity and, as a result, an almost two-order-of-magnitude increase of the statistics as compared with the previous-generation experiment Belle.

V. THE POLE POSITIONS OF THE
Z b , Z 0 b AND W bJ STATES A. Extracting the poles in a multichannel scattering problem In this section we employ the approach developed above to predict in a parameter-free way the pole positions for the spin partners W bJ (J ¼ 0, 1, 2) with the quantum numbers J þþ . For the pole search in the complex energy plane we follow the approach of Refs. [15,20] and stick to the foursheet Riemann surface corresponding to two elastic channels-either the BB and B ÃBÃ channels in the case of J PC ¼ 0 þþ or the BB Ã and B ÃBÃ ones for all other quantum numbers. All inelastic thresholds are remote and their impact on the poles of interest, which are located near the elastic thresholds, is minor-see the discussion in Ref. [15]. Then, for two coupled channels with the thresholds split by the mass difference Δ, the four-sheeted Riemann surface can be mapped onto a single-sheeted plane of a new variable, which is traditionally denoted as ω [48,49], via the relations 1 Then, the energy defined relative to the lowest threshold of the two, where μ 1 and μ 2 are the reduced masses in the first and second elastic channels labeled as (1) and (2), respectively. Specifically, in the 0 þþ channel, while in the channels 1 þþ and 2 þþ one has Then the one-to-one correspondence between the four Riemann sheets in the E-plane (denoted as RS-N, where N ¼ I, II, III, IV) and various regions in the ω-plane read RS-IðþþÞ∶ Im k 1 > 0; Im k 2 > 0; RS-IIð−þÞ∶ Im k 1 < 0; Im k 2 > 0; RS-IIIð−−Þ∶ Im k 1 < 0; Im k 2 < 0; RS-IVðþ−Þ∶ Im k 1 > 0; Im k 2 < 0; where the signs in the parentheses correspond to the signs of the imaginary parts of the momenta k 1 and k 2 . These regions in the ω-plane are depicted in Fig. 11. The thick solid line corresponds to the real values of the energy lying on (physical) RS-I. It is easy to see that the physical region between the two thresholds corresponds to jωj ¼ 1, with both ReðωÞ and ImðωÞ positive, and the thresholds at E ¼ 0 and E ¼ Δ are mapped to the points ω ¼ AEi and ω ¼ AE1, respectively. The nomenclature of the Riemann sheets as described above is relevant for a two-channel situation while in the presence of additional channels it needs to be generalized. As an illustration, consider a three-channel case with two elastic and one inelastic channel. In this case, the threechannel complex omega plane can be schematically viewed as a two-sheeted ω-plane with the sheets connected by analyticity through the inelastic cut (see Fig. 12) where, in line with the two-channel case above, the sheets are labeled TABLE IV. The ratios of the individual widths for the elastic and inelastic channels in the ϒð10860Þ radiative decays via the W bJ states relative to the sum of all individual partial widths for a given J (all ratios in each line add up to unity) obtained for the pionful fit 2 (fit G in Ref. [15]). poles residing in the right ω-plane do not meet this criterion and, therefore, can be safely disregarded. For example, in order to reach I up ðþþþÞ starting from the lower domain of the RS-I [I low ðþþþÞ in the right ω-plane-see Fig. 12] one would need to travel a long way around the inelastic branch point. Among the poles residing in the left ω-plane those on the sheets ð−−þÞ and ð−−−Þ are the most important ones since they are the closest to the physical region of the real energy (the fat black line in Fig. 12). A pole on the sheet ð−−þÞ located in the vicinity of the lower elastic threshold (the point ω ¼ i) will result in a significant near-threshold distortion of the line shapes while if the pole is located deeper in the complex plane (but still on the same sheet) it will show up as a clear resonance peak. The poles residing on the sheet ð−−−Þ will manifest themselves in the observables in exactly the same manner but with respect to the upper threshold. Next in importance are the poles residing on the sheets ðþ−þÞ and ðþþ−Þ. They cannot generate resonance humps in the line shapes above threshold but can significantly enhance the cusplike structure at threshold provided these poles reside not too deep in the complex plane. It is also important to notice that imposing constraints from unitarity and analyticity on the scattering amplitude T requires that This implies that if there is a pole at ω ¼ ω 0 in the left (main) omega plane there must be also a pole at −ω Ã 0 in the right omega plane-see Fig. 12. This mirror pole, however, has no physical significance, as already explained.
Generalization of the logic discussed above to a larger number of channels is straightforward if one bears in mind that, as before, only one ω-plane sheet containing the domain of the physical energies on RS-I is relevant and that the signs of the imaginary parts for all remote inelastic channels coincide.
For the case at hand, we arrive effectively at a threechannel (two elastic plus an effective inelastic) problem for the quantum numbers 1 þ− , 1 þþ and 0 þþ while for 2 þþ , where all three elastic channels are present, the effective problem contains four channels.  11. The unitary-cut-free complex ω-plane for the two elastic channels obtained from the four-Riemann-sheeted complex energy plane by the conformal transformation (75). The eight regions separated by the unit circle and by the two axes correspond to the upper and lower half-planes (see the subscripts "up" and "low") in the four Riemann sheets of the energy plane denoted as RS-N with N ¼ I, II, III, IV [48,50]. The bold line indicates the physical region of a real energy E [48].

B. Poles of the Z b 's and their spin partners
In the vicinity of a pole located at M ¼ M R α the elastic scattering amplitude T αα ðM; p; p 0 Þ given in Eq. (62) takes the form where the energy M is defined in Eq. (66) and g 2 α and M R α stand for the residue and the pole position in the channel α, respectively.
The most relevant poles for the case of the quantum numbers J PC ¼ 1 þ− are collected in Tables V-VII, together with the residues at these poles. As was explained in detail above, we regard a pole as relevant if it has a short (compared with the splitting between the nearest elastic thresholds δ) path to the physical RS-I and as such affects the form of the line shapes. In the pionful fits, the poles representing the Z b ð10610Þ and Z b ð10650Þ states inhabit the sheets ð−−þÞ and ð−−−Þ, respectively, and, according to the logic discussed above, reveal themselves in the line shapes as peaks above thresholds.
Also, in Tables V-VII, we present the relevant poles (counted relative to the nearby elastic thresholds) and the corresponding residues predicted for the spin partners in the 0 þþ , 1 þþ and 2 þþ channels. For the pionful fit 2 regarded here as the most reliable calculation, the poles in all channels reside on the sheets closest to the upper domain of the physical RS-I. The poles for the Z b and W 0 b0 are located just in the vicinity of the BB Ã and B ÃBÃ threshold, respectively, and show up as near-threshold distortions in the line shapes (see the black solid lines in Figs. 5 and 7). Meanwhile, the poles representing the other states are shifted from the respective thresholds to the complex plane by about 10-15 MeV and manifest themselves as humps (Z 0 b and W b0 ) or pronounced above-threshold peaks (W b1 and W b2 ) (see the black solid lines in Figs. 7-9). The shift of the pole positions in the pionful fit 2 as compared with the contact fit, where all poles correspond to virtual states, appears mainly due to the pion dynamics-this effect is fully in line with the findings of Ref. [14].

C. Uncertainty estimate
Uncertainties of the poles and residues given in Tables V-VII correspond to a 1σ deviation in the parameters of the fit from the central values shown in Table III. The source of this uncertainty is the experimental errors in the data. The ω-plane sheet closest to the physical region of a real energy indicated by the black bold line. In this ω-plane, the inelastic momentum for the particles with the masses m in1 and m in2 is related to the energy as The ωplane sheet distant from the physical region. Here the inelastic momentum is related to the energy as Transitions from one omega plane to the other are possible through the right-hand cut from the inelastic channels denoted by the fat pink line in both panels.
As for the theoretical uncertainty, it can be estimated as the maximum of the two errors from the truncation of the EFT expansion at a given order (for the discussion of the truncation errors in the NN sector see, e.g., Ref. [51]) and from the cutoff variation. To explain the truncation error method, let us introduce an observable quantity X ðνÞ ðQÞ calculated to a given order ν in the EFT expansion in the momentum Q, where Q ∼ p typ ¼ 0.5 GeV and Λ h ∼ 4πf π ≃ 1 GeV with f π denoting the pion decay constant; fα n g are the expansion coefficients with α 1 ¼ 0 since there are no operators at the order Q. Then, assuming that the expansion (80) converges, the error at the given order ν is expected to come from the first neglected chiral order; that is, it should scale as χ νþ1 unless the coefficient α νþ1 vanishes. In the latter case, the uncertainty is estimated based on the nonvanishing result at the order χ νþ2 , and so on. For example, the observable at LO (ν ¼ 0) and its truncation error read where we used that α 1 ¼ 0.
Although the poles are not observed directly, they manifest themselves in observable quantities such as line shapes, and thus the truncation error method is expected to work for them too. To estimate the truncation error at LO (ν ¼ 0), we compare the results calculated explicitly at the orders ν ¼ 0 (pionful fit 1) and ν ¼ 2 (pionful fit 2) to find that the truncation error for the W b2 does not exceed 5 MeV while it is only about 1 MeV for the near-threshold state W 0 b0 as well as for the Z b and Z 0 b states. On the other hand, the truncation error at LO for the states W b1 and W b0 is of the order of 15 MeV which does not look unnatural either given the large expansion parameter of the pionful EFT. It also needs to be emphasized that the chiral expansion for the W b0 state might converge slower than expected from our power counting, because this state lies by δ ¼ 45 MeV lower than the energy region used in the fits for the Z b 's. It remains to be seen how the pole position for this state is affected by the inclusion of higher-order interactions.
The truncation error method becomes particularly useful when the results at least at several chiral orders are calculated explicitly, which is not yet feasible. Still, from the results presented above one can conclude that the partner states W b2 and W 0 b0 , both residing near the B ÃBÃ threshold, indicate a very good stability with respect to the inclusion of higher-order interactions. Since the truncation error estimate for the NLO results (pionful fit 2) is not possible at present (it would require a complete N 2 LO calculation) we rely on naturalness to provide a rough TABLE V. The pole positions and the residues g 2 [see the definition in Eq. (79)] in various S-wave B ðÃÞBðÃÞ channels for the contact fit. The energy E pole is given relative to the nearest open-bottom threshold quoted in the third column, so that it is one of the energies E n (n ¼ 1, 2, 3) defined in Eq. (66). The Riemann sheet (RS) is defined by the signs of the imaginary parts of the corresponding momenta (quoted in the columns 4-7); a missing sign indicates that this channel is uncoupled. Uncertainties correspond to a 1σ deviation in the parameters allowed by the fit to the data in the channels with J PC ¼ 1 þ− where the Z ð0Þ b states reside [15]. For the estimate of the theoretical uncertainties see Sec. V C. The poles are calculated for the cutoff Λ ¼ 1 GeV.

J PC
State Threshold Im p in Im p BB Im p BB Ã Im p B ÃBÃ E pole w.r.t. threshold (MeV) Residue at E pole This estimate gives roughly a twice larger uncertainty for the pionful fit 2 than the cutoff variation.

VI. SUMMARY AND CONCLUSION
In this paper we address the properties of the spin partners W bJ with the quantum numbers J þþ (J ¼ 0, 1, 2) of the bottomoniumlike states Z b ð10610Þ and Z b ð10650Þ. We employ the EFT approach developed previously in Ref. [15] and fix all unknown low-energy constants and couplings from the data on the line shapes in the elastic and inelastic channels for the negative C-parity states Z b and Z 0 b . After that, the same EFT approach consistent with requirements from unitarity, analyticity and HQSS is employed to predict in a parameter-free way the line shapes of the positive C-parity spin partner states W bJ in the corresponding elastic [B ðÃÞBðÃÞ ] and inelastic [η b ðnSÞπ and χ bJ ðmPÞπ] channels.
Because of the positive C-parity the W bJ 's should be produced in the radiative decays of the vector bottomonium ϒð10860Þ. It is argued that the production operator which involves the tree-level and one-loop contributions behaves as a smooth function of the energy in the near-threshold region of interest here. Therefore, in agreement with the Watson's theorem, the energy dependence of the line shapes can be predicted based on the strong interaction amplitudes between the heavy mesons in the final state. These amplitudes contain the poles in the vicinity of the thresholds which are associated with the excitation of the W bJ states. In addition, the ratios of the partial branchings to all aforementioned elastic and inelastic decay channels of the W bJ partners come as predictions of our approach, since the overall normalization constant from the production operator drops out in these ratios in the HQSS limit.
With a multichannel amplitude at hand which possesses the correct analytic structure we extract the poles of the amplitude in the complex energy plane and its residues at these poles for all four partner states and find that our most advanced pionful analysis of the data on the Z b 's (the pionful fit 2) is consistent with all W bJ 's being abovethreshold resonances. In contrast to this, in the pionless approach, all W bJ 's appear as virtual below-threshold states. Since these two scenarios reveal themselves differently in the line shapes (cf. the threshold cusp versus the above-threshold hump in the inelastic channels in Fig. 10), the experimental data in various channels relevant for the W bJ 's should provide key information about the role of the pion dynamics for the system at hand.
The uncertainties in the pole positions are estimated and discussed in detail. The errors accounted for in this work come from (i) a statistical 1σ deviation in the parameters allowed by the fit to the data, (ii) truncation of the EFT expansion at a given order and (iii) the cutoff variation. Although the evaluated uncertainties appear to be of a natural size a better estimate of the truncation error would be very desirable. That would call for the inclusion of the two-pion exchange contributions to the elastic potentials at next-to-leading order which does not involve any new parameters.
It is important to notice that the developed chiral EFT approach relies on the most general effective Lagrangian and, within the domain of its validity (the energy region of the order of δ ¼ 45 MeV around the BB Ã threshold), can be used to make model-independent predictions for nearthreshold exotic states (line shapes, poles and residues) regardless of particular origins of the short-range interactions involved. The importance of bare poles (for example, associated with the compact tetraquarks) in the regime of applicability of the EFT would bring additional scales into the problem and result in a strong energy dependence of the short-range interactions in the EFT. This would spoil the convergence of the EFT derived above, which is, however, not supported by the results of our studies. Indeed, the analysis performed in Ref. [15] and also here shows that the derivative expansion for the contact terms is consistent with the existing experimental data for the Z b states, which are, therefore, well described without the necessity of including compact tetraquarks explicitly.
Once the parameters of a resonance are extracted, its nature can be inferred from the analysis of the pole positions and the residues which provide the coupling of a state to the continuum channels. In particular, we observe that, as expected, for each Z b or W bJ state only one nearthreshold pole in the complex momentum plane defines its properties (see Table VII for the corresponding pole in the energy plane). According to the pole counting rules suggested in Ref. [52], this indicates dynamical origins of the studied resonances-there would have been two nearly symmetric poles for a compact quark state. This conclusion is further supported by the large dimensionless (of the order of unity) coupling constants extracted from the residues at the poles for both Z b 's and their spin partners which result in the large partial branchings of their decays to the open-flavor channels-see Table IV. We believe that all these manifestations together provide a strong argument in favor of the dominating molecular component in the wave functions of these states. We conclude by stating that, although the electromagnetic fine structure penalty suppresses the probability of the W bJ 's production in radiative decays of the ϒð10860Þ by two orders of magnitude compared with the Z ð0Þ b production in the one-pion decays of the ϒð10860Þ, this suppression is to be overcome by the large statistics anticipated for the Belle-II B-factory. We, therefore, expect the spin partners of the Z b states to be copiously produced in this experiment.
where the contact terms proportional to the LECs D 10 and D 11 contribute to S-wave interactions while the term D 12 gives rise to the S-D transitions. As explained in Ref. [15], we are only interested in the S-S and S-D transitions for the B ðÃÞBðÃÞ scattering, so that all terms of the kind ∝ ∇ i H † ∇ j H contributing to P waves are dropped.
Further, we define the combinations