The line shapes of the $Z_b(10610)$ and $Z_b(10650)$ in the elastic and inelastic channels revisited

The most recent experimental data for all measured production and decay channels of the bottomonium-like states $Z_b(10610)$ and $Z_b(10650)$ are analysed simultaneously using solutions of the Lippmann-Schwinger equations which respect constraints from unitarity and analyticity. The interaction potential in the open-bottom channels $B^{(*)}\bar{B}^{*}+\mbox{c.c.}$ contains short-range interactions as well as one-pion exchange. It is found that the long-range interaction does not affect the line shapes as long as only $S$ waves are considered. Meanwhile, the line shapes can be visibly modified once $D$ waves, mediated by the strong tensor forces from the pion exchange potentials, are included. However, in the fit they get balanced largely by a momentum dependent contact term that appears to be needed also to render the results for the line shapes independent of the cut-off. The resulting line shapes are found to be insensitive to various higher-order interactions included to verify stability of the results. Both $Z_b$ states are found to be described by the poles located on the unphysical Riemann sheets in the vicinity of the corresponding thresholds. In particular, the $Z_b(10610)$ state is associated with a virtual state residing just below the $B\bar{B}^{*}/\bar B{B}^{*}$ threshold while the $Z_b(10650)$ state most likely is a shallow state located just above the $B^*\bar{B}^{*}$ threshold.


Introduction
In the last decade, numerous new hadrons have been observed in the charmonium and bottomonium energy region. Of special interest are those that cannot be accommodated by the simple quark-model picture and which are, therefore, referred to as exotic hadrons. For example, the charged states Z ± b (10610), Z ± b (10650) [1], Z ± c (3900) [2,3], Z ± c (4020) [4], Z ± (4430) [5][6][7][8] which, amongst other channels, decay into quarkonium states and a pion and cannot be conventionalQQ (with Q denoting the heavy quark) mesons as their minimal quark contents is four-quark.
At present, both a tetraquark structure [12][13][14] and a hadronic molecule interpretation [15][16][17][18][19][20][21][22][23] are claimed to be consistent with the data for these two exotic states. Their proximity to the BB * and B * B * thresholds together with the fact that those are by far the most dominant decay channels of the Z ± b (10610) and Z ± b (10650), respectively, provides a strong support for their molecular interpretation. While some works try to explain the Z b 's as a simple kinematical cusps [24], it was demonstrated in Ref. [25] that the narrow structures in the elastic channels necessitate near-threshold poles. The general argument presented there was supported by the explicit analyses presented in Refs. [26,27]. Also the analysis presented in this work finds poles near the thresholds in conflict with the claims of Ref. [24].
The literature on the hadronic molecule scenario for the Z b states is already very rich: hadronic and radiative decays are studied in Refs. [28][29][30][31][32][33][34], the contribution of the two Z b states to other processes are considered in Refs. [35][36][37], the heavy-quark spin partners of the Z b 's are discussed in Refs. [21,[38][39][40][41][42][43][44], the line shapes and pole positions are addressed in Refs. [16,26,27,45,46], and the problem of three-body universality in the B-meson sector is investigated in Ref. [47]. For a recent review of the theory of hadronic molecules we refer to Ref. [48]. Since both the Z ± b (10610) and Z ± b (10650) contain a heavy bb pair, it is commonly accepted that the constraints from the heavy-quark spin symmetry (HQSS) should be very accurate for these systems. For instance, HQSS allows one to explain naturally the interference pattern in the inelastic channels Υ(nS)π and h b (mP )π [15]. On the contrary, there is a long-lasting debate in the literature about the role of the one-pion exchange (OPE) interaction played in the formation of hadronic molecules. In a pioneering work [49] a vector meson exchange was proposed as a key ingredient of the potential between a heavy meson and a heavy antimeson. In contrast to this the OPE was used in analogy to the deuteron in Refs. [50,51] to predict the existence of the X(3872). The model was further elaborated and extended to other channels in Refs. [52,53], but it was criticised in Refs. [54][55][56], where amongst other issues the potential importance of the three-body dynamics was stressed.
A more advanced approach for the X(3872) and other molecular candidates is based on an effective field theory (EFT) treatment which incorporates both short-range interactions parameterised by a contact term and long-range interactions due to the OPE. 2 There exist two competing EFT approaches: one that treats the pions as a perturbation on top of the nonperturbative contact interaction -the so-called X Effective Field Theory (X-EFT), see, for example, Refs. [58,59], and the other in which both the contact and the OPE interactions are iterated to all orders fully nonperturbatively [60]. Both approaches can be generalised naturally to the b-quark sector. In particular, based on the relatively weak coupling of the pion to heavy fields in Refs. [39,59] a power counting scheme is proposed within the X-EFT framework to conclude that the central OPE can be included perturbatively in the heavy-quark sector. On the other hand, it was noticed in Refs. [44,61] that the most prominent contribution from the OPE stems not from the S-S but from the S-D coupled-channel transitions generated by tensor forces and that reliable predictions for the spin partners of the X and Z b 's can only be made if the OPE interaction is included nonperturbatively. 3 For example, due to the OPE both the mass and the width of the 2 ++ partner of the X(3872) experience a considerable shift as compared to the X-EFT based predictions made in Ref. [63] and the properties of the spin partners of the Z b 's are also quite sensitive to the pion exchange [44] although to a lesser extent than in the c-quark sector.
Within recent phenomenological studies, it was claimed in Ref. [64] that, near the Swave open-bottom thresholds BB * andB * B * , the effect of the OPE on the line shapes of the Z b 's can be as large as 30%. On the other hand, it was advocated in Ref. [65] that the OPE gets cancelled by the one-eta exchange (OEE), so that in practice it is irrelevant for the formation of the molecular states. Thus, we take all those controversies as a motivation to investigate in detail the role of the OPE for the properties of the Z b states. In particular, in this work we concentrate on the line shapes and pole locations.
The fact that the full OPE has to be incorporated into a coupled-channel approach to the B ( * )B( * ) system should be obvious already from the momentum scales involved: For energies near the B * B * threshold, the on-shell relative momentum in the B * B channel is as large as MeV denotes the B * -B mass difference, with m B and m B * being the B and B * meson mass, respectively. While the D waves in the OPE are indeed suppressed for momenta much smaller than the pion mass (p m π ), in the opposite regime of relevance here (p typ m π ), S-D transitions mediated by the OPE turn out to be as important as S-S transitions. The appearance of such a scale as p typ implies that the convergence of the EFT approach is controlled by the expansion parameter where Λ h denotes a typical hadronic scale of the order of 1 GeV. As a consequence, the EFT might converge relatively slowly and higher-order interactions such as O(p 2 ) contact terms and HQSS violating contact interactions might play a role for the line shapes too. In this paper, we investigate the impact of the OPE, OEE and the higher-order interactions on the line shapes for the near-threshold states Z b (10610) and Z b (10650). In order to appropriately account for the scale p typ , we iterate the resulting potential to all orders. In Refs. [26,27], the relevant data were analysed using a parameterisation for the line shapes in both elastic and inelastic channels consistent with unitarity and analyticity based on the leading (O(1)) S-wave short-range interactions only. In this paper, we reanalyse the same set of data using a direct numerical solution of the Lippmann-Schwinger Equations (LSEs). This allows us not only to check the validity of the approximations made in Refs. [26,27] in order to derive self-consistent closed-form analytic expressions but also to include nonperturbatively the pion exchange and other contributions from the scattering potential.
The paper is organised as follows. In Sect. 2 we discuss the effective potentials in the system at hand. In particular, in Subsect. 2.1, we start from the theory with purely contact transition potentials between the elastic and inelastic channels. Then, in Subsect. 2.2, we discuss the one-pion and one-eta exchange interactions between the B ( * ) mesons and include them on top of the contact interactions. The resulting Lippmann-Schwinger equations are derived and solved numerically in Sect. 3. Different fitting strategies and the results of the data analysis are presented in Sect. 4. Sect. 5 is devoted to searches of the poles in the complex plane which describe the Z b states. We summarise in Sect. 6.

Effective potentials
The effective potentials between two heavy mesons, which enter LSEs, contain local contact terms and the contributions from the lightest pseudoscalar Goldstone boson octet. We start from a discussion of the short-range contributions parameterised by the contact potentials without derivatives and with two derivatives. Then we discuss in detail the one-pion and one-eta exchange potentials relevant for the study. In what follows, we stick to the labels introduced previously in Refs. [26,27]: The N e open-flavour channels (qb)(bq) (here q is a light quark) are denoted by greek letters α, β and the N in hidden-flavour channels (bb)(qq) -by latin letters i, j. Elementary poles which would represent compact quark compounds are not considered because the minimal quark contents of the isovector Z b states is fourquark. In the system at hand, the elastic channels are BB * and B * B * and the inelastic channels are Υ(nS)π and h b (mP )π with n = 1, 2, 3 and m = 1, 2. Therefore, N e = 2 and N in = 5.

Short-range contributions
Following the notation introduced above, the full pionless potential takes the form of a (2N e + N in ) × (2N e + N in ) matrix, where the basis vectors are denoted as with the letters S or D in square brackets standing for the orbital angular momentum L in the corresponding elastic channel. Indeed, the BB * and B * B * systems with the quantum numbers J P C = 1 +− can be either in the 3 S 1 or in the 3 D 1 state.
In what follows we assume that all inelastic channels only couple to the S-wave elastic ones as their couplings to the D-wave elastic channels are suppressed by a factor p 2 typ /m 2 B 1, since the transitions are driven by the exchange of a B meson -for details we refer to App. B. Thus, we set v iα = 0 for α = BB * [D] and B * B * [D]. Furthermore, the transition potentials between the α-th elastic channel in an S wave and the i-th inelastic channel can be parameterised through the coupling constants g iα , where k i and l i are the momentum and the angular momentum in the i-th inelastic channel, respectively. The nonvanishing (S-wave) couplings g iα are subject to the HQSS constraints [26,27,66], where, as before, n = 1, 2, 3 and m = 1, 2. Therefore, in what follows, as long as we discuss the results in the HQSS limit, only the coupling constants for the BB * channel are quoted in the form Furthermore, following the arguments from Refs. [26,27], we neglect the direct interactions in the inelastic channels, since their thresholds are located far away from the relevant energy region and the interaction of light states withQQ states is suppressed, thus setting v ji (k , k) = 0 for all i's and j's -see Eqs. (2.1) and (2.2) -that allows us to disentangle the inelastic channels from the elastic ones and to reduce the effect of the former to just an additional term in the effective elastic-to-elastic contact transition potential, with m H i (m h i ) denoting the mass of the heavy(light) meson in the i-th inelastic channel and M being the total energy of the system. Furthermore, the inelastic momentum is defined as where λ(m 2 1 , m 2 2 , m 2 3 ) is the standard triangle function. At leading order O(p 0 ) (LO), the short-range elastic potential v αβ in the strict HQSS limit consists of two momentum-independent contact interactions [15,28,38,41,67] (see also Appendix A for further details), so that for the basis vectors α (see Eq. (2.3)) it can be written as In what follows, we will also investigate the influence on the line shapes of HQSS violating corrections and next-to-leading order (NLO) contact interactions with two derivatives. To this end, we extend the elastic potential by including three additional contact terms proportional to D d , D l and D SD at the order O(p 2 ) (see Appendix A for the corresponding Lagrangian densities), where the first two low-energy constants (LECs) give rise to the S-S transitions while D SD contributes to the S-D (and D-S) ones. In addition, we introduce the leading HQSS violating contact interaction which contributes to the S-S diagonal transitions, so that for the resulting elastic potential we arrive at (2.10)

One-pion and one-eta exchange potentials
The potentials for the OPE and the OEE at LO can be derived from the effective Lagrangian, Tr σ · u abHbH † a + h.c., (2.11) where the superfield H a , describes the heavy-lightqQ mesons (its transformation properties are discussed in Appendix A), the axial current is u = −∇Φ/f π + O(Φ 3 ), the matrix for the pseudoscalars reads (here only the matrix responsible for the SU(2) subspace relevant here is retained for simplicity) 13) and the pion decay constant is f π = 92.4 MeV [68]. Then, the Lagrangian of the Φ-field coupling to a pair of heavy-light mesons reads, to leading order, (2.14) In the strict heavy-quark limit, the dimensionless coupling g Q does not depend on the heavy quark flavour Q, so that g b = g c . The latter can be extracted from the D * partial decay width, where k is the three-momentum in the final state measured in the rest frame of the decaying particle and m D and m D * are the D-and D * -meson mass, respectively. Then, from the width Γ(D * + → D 0 π + ) = 56.46 keV [68] one extracts which agrees within 10% with the recent lattice QCD determination of the B * Bπ coupling constant [69] and which will be used in all calculations below.
In this section only the OPE potential will be discussed in detail since the OEE potential can be obtained straightforwardly from the expressions for the OPE by replacing (i) the flavour factor +1 (combining the isospin coefficient −1 for the OPE in the isovector channel and the negative C-parity factor −1) by −1 for the OEE; (ii) the pion mass, m π , by the η mass, m η ; (iii) the pion coupling constant g b by the η coupling constant g b / √ 3 -see, for example, Ref. [70] and Eq. (2.13) above. Further details can be found in Appendix C.
In the framework of the Time-Ordered Perturbation Theory (TOPT), the OPE potential acquires two contributions depicted in Fig. 1 with M being the total energy of the system and E x being the energy of the particle x.
The incoming and outgoing momenta are denoted as p and p , and the Cartesian indices n and n contract with either the B * (B * ) polarisation vectors or with the vector product of polarisation vectors, B * ×B * . In the presence of pions, no additional coupled-channel transitions but those discussed in Sect. 2.1 are possible for the B ( * )B * systems with the quantum numbers 1 +− . Therefore, using the elastic basis vectors introduced in Sect. 2.1 -see Eq. (2.3) -the OPE potential in this channel can be written as 4 where, in each matrix element V π λλ LL (M, p, p ), the index λ(λ ) = 1, 2 labels the particle channel (BB * = 1, B * B * = 2) and L(L ) stands for the angular momentum of the state λ(λ ). Then, The functions ∆ π λλ k (k = 0, 1, 2, 3) are defined as with the contributions of the two TOPT orderings -see Fig. 1 -given by where p π = p 2 + p 2 − 2pp x, and the mass matrices for the intermediate particles read Note that, in the static limit, that is, when the recoil corrections are neglected, the functions V π λλ 1 , V π λλ 2 , and D π λλ a (p, p , x), D π λλ b (p, p , x) do not depend on the particle-channel indices (λ and λ ), so that in this approximation the dependence on the latter comes entirely from the coefficients in Eq. (2.21).

Numerical solutions of the coupled-channel Lippmann-Schwinger Equations
With the OPE and OEE interaction included and decomposed in partial waves, the full interaction potential in the elastic channels reads where the effects of the inelastic channels are contained in V CT αβ (M, p, p ) as given by Eq. (2.7) and the channel indices are defined by Eqs. (2.2) and (2.3).
We work in terms of the production amplitudes rather than the scattering amplitudes which are more convenient given that we aim at fits for the invariant mass distributions measured in the Υ(10860) decays. Thus, the system of the LSEs reads Here F α (M, p) and U α (M, p) denote the Born and the physical production amplitude of the α-th elastic channel from a point-like source, respectively. We include those source terms in S waves only. Note that a D wave contribution to the source could be re-expressed as an additional energy dependence of the source, which, however, is not required by the data. This is also in line with the results of Ref. [72] from which one concludes that the two-step process Υ(10860) → B * B ( * ) → B * B ( * ) π with the pion emitted from the B-meson line, that provides the energy dependence, is suppressed as compared to a point-like source. The Green's function for a two-heavy-meson intermediate state reads where m α th stands for the α-th elastic threshold and µ α is the reduced mass in the channel α. Other components of the multichannel amplitude responsible for the 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 [26,27]. In particular, for the i-th inelastic channel in the final state we have where λ(m 2 1 , m 2 2 , m 2 3 ) is the standard triangle function. Expression (3.4) is based on the assumption that the data are dominated by the Z b (10610) and Z b (10650) poles which emerge from B ( * )B * dynamics, so that the Born amplitudes F i (M, p) coming from the inelastic sources can be safely neglected. While this is well justified for the h b (mP )π channels, since the Z b poles are necessary for the change in the heavy quark spin, it appears to be unjustified for the heavy-spin-conserving Υ(nS)π channels. However, to fully control the interplay of the source term and the resonance terms one would need to properly include the ππ interaction which goes beyond the scope of this work. 5 This is why, in what follows, we do not include the line shapes in the Υ(nS)π channels into the fit.
In terms of the production amplitudes the expressions for the differential widths in the elastic and inelastic channels read (see Refs. [26,27] for the derivation) 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 (i-th inelastic) channel in the rest frame of the B * B( * ) (Υ(nS)π/h b (mP )π) system. Then, the total branching in an elastic or inelastic channel x is defined as where For simplicity, we define the branching fractions (BFs) relative to the BB * π channel, The χ 2 function to be minimised in the fitting process is built as where y exp 's are the experimental distributions given in the form of histograms (the sums in n's run over bins), δ's denote the errors, and the normalisation factors N x are additional (auxiliary) fitting parameters since the data are presented in arbitrary units. Both line shapes and total branchings are used in the fit for the elastic B ( * )B * and the inelastic h b π channels while for the inelastic channels Υπ only the total branchings can be used in the one-dimensional analysis.
Since the production proceeds via pion emission from the Υ(10860) bottomonium, the corresponding Born amplitudes are subject to a HQSS constraint similar to the one given in Eq. (2.5) for the coupling constants, that is, [26,27]. Then, without loss of generality, one can set to fix the overall normalisation of the amplitudes which in any case drops out from the branchings and can be absorbed by the unimportant factors N in the χ 2 function -see Eqs. (3.7) and (3.10).
We may now perform the angular integrations in Eq. (3.2) to reduce the three-dimensional integral equation to a one-dimensional equation, To render the integrals well defined we introduce a sharp ultraviolet cut-off Λ which needs to be larger than all typical three-momenta related to the coupled-channel dynamics. Unless stated otherwise, for the results presented below we choose Λ = 1 GeV. The cut-off dependence of our results will be discussed in Sect. 4.2. When the energy goes above the threshold of the intermediate channel β, the on-shell three-momentum p β is real allowing for a singular integrand at q 2 = p 2 β . In this case one subtraction at the point q = p β is implemented to stabilise the numerical result, Since this equation is valid in the whole complex energy plane, it is also used to find the resonance poles in Sect. 5. The masses of particles used in the calculations are [68] m B = 5279 MeV,

Fit schemes and results
In this section we fit the line shapes of the two Z b states in both elastic (B ( * )B * ) and inelastic (h b (mP )π, m = 1, 2) channels. As was already mentioned above, the line shapes in the inelastic Υ(nS)π (n = 1, 2, 3) channels cannot be included into the fit yet, since the BF, % Exp. 50 ± 10 0.6 ± 0.  Table 1. The branching fractions (BFs), in per cent, for the elastic and inelastic channels in the Υ(10860) decays via the two Z b states relative to the BB * π channel for which the BF is set to unity [1,11,73]. The branching fractions within fit schemes A and G are also listed as a comparison.
data contain a significant contribution driven by the two-pion final state interaction that we cannot include straightforwardly in the present approach. In addition, the analysis has to be multidimensional. Meanwhile, the partial branchings for all measured elastic and inelastic channels are included in our analysis -see Eq. (3.10) for the formula for χ 2 and Table 1 where the branchings are quoted relative to the BB * π channel. In the branching fractions the interference terms with crossed channels are effectively included via scaling factors gauged in test calculations. It turns out that those scaling factors are very close to 1 in the h b π channels, but are 5%, 10%, and 20% for the partial BF's in the Υ(1S)π, Υ(2S)π, and Υ(3S)π channels, respectively.

The role of pion dynamics, HQSS violation and higher-order interactions
In this section, we investigate the impact of various effects, included in this work for the first time, on the line shapes. To gain a proper insight on the role of these effects, we include them one by one. We, therefore, consider eight different schemes: • Scheme A: Only the O(p 0 ) S-wave contact potentials are considered.
• Scheme B: As Scheme A with OPE added, however, only in S waves.
• Scheme C: As Scheme B, but with D-wave OPE included.
• Scheme D: As Scheme C but allowing for a sizeable HQSS violation.
• Scheme E: As Scheme C, but with the O(p 2 ) D SD contact potential.
• Scheme F: As Scheme C, but with the O(p 2 ) D d , D f , D SD contact potentials.
• Scheme G: As Scheme F, but with OEE included in addition.
• Scheme H: As Scheme A, but with the O(p 2 ) D d , D f , D SD contact potentials in addition.
The parameters of the fits for the eight schemes described above are listed in Tables 2  and 3. In what follows, we discuss the quality of the fits and draw conclusions on the role played by various effects. Scheme The fit results for Schemes A, B and C are shown in Fig. 3 by the solid black, red dotted and blue dashed curves, respectively. The line shapes of Scheme A are basically identical to those based on the parameterisation proposed in Refs. [26,27,66]. In particular, also here the diagonal and off-diagonal matrix elements of the direct interaction potential satisfy the strong inequality |C f | |C d | which ensures that the transitions between the two elastic channels are suppressed. It has to be noticed, however, that the parameters in Scheme A of the present paper cannot be directly compared with those from Ref. [27] because of a different regularisation.
Scheme B, with the S-wave dynamic OPE included, provides a fit of the same quality as the one in Scheme A. The resulting line shapes are quite similar for both fits which is consistent with the claim made in Ref. [27] on a moderate role played by the OPE. On the  contrary, the claim of Ref. [64] that already the S-wave OPE changes the line shape up to 30% near threshold is not supported by our fit B. It is instructive to identify the most relevant differences between the approaches employed in Ref. [64] and in this work: First, the off-diagonal terms connecting the BB * with the B * B * channel and neglected in Ref. [64] turn out to be large in the OPE potential and even exceed the diagonal terms by roughly a factor of 2 -see Eq. (2.21). Secondly, the argument of Ref. [64] is based on a perturbative treatment of the OPE while in this work the pions are treated nonperturbatively and all parameters are re-fitted to the data. Indeed, before drawing conclusions on the role of longrange interactions the short-range part of the OPE, intimately connected to the contact interactions [57], needs to be appropriately renormalised which is achieved in this work by refitting the contact terms to the line shapes. After the re-fit we find that the central S-wave OPE can be absorbed to a very large extent into a re-definition of the contact interactions which is in line with the observation made in Ref. [44].
The fit for Scheme C demonstrates that D waves in the OPE can play a non-negligible  role for the line shapes -this is most clearly visible in the BB * spectra in Fig. 3 where now a bump appears around the B * B * threshold. This result should not come as a surprise given the large momentum scale p typ , defined in Eq. (1.1), introduced by the splitting between the elastic thresholds, which enhances the D waves and, in particular, the contribution from the S-D transitions. These findings are in line with the claims made in Refs. [44,61]. Although the fit for Scheme D is essentially consistent with the B * B * distribution as well as with the distributions in the inelastic channels, the shape of the BB * spectrum distorted by the D-wave OPE is not supported by the data. The observation that the experimental line shapes do not exhibit a hump structure around the B * B * threshold was related in Ref. [74] with a possible existence of the light-quark symmetry in QCD. In what follows, we investigate two different variations of the potential aiming at an improved description of the data.
In Fig. 4 we demonstrate the impact of various higher-order interactions on the line shapes. In particular, in the fit for Scheme D (green dot-dashed curves) we release the HQSS breaking parameters in the elastic and inelastic potentials as well as in the production vertex (that is, ξ Υ(nS) , ξ h b (mP ) and ) and allow them to deviate up to 50% from the HQSS predicted values. This gives However, in spite of a significant HQSS breaking allowed in the fit, the resulting distribution does not show a qualitative improvement leaving, in particular, the bump structure around the B * B * threshold nearly unchanged.
On the other hand, the inclusion of a single O(p 2 ) contact interaction D SD between the S-wave and D-wave elastic channels (Scheme E) improves the fit considerably (see the black dotted curves) yielding χ 2 /d.o.f = 1.11, as given in Table 2. As required by the data, this term is fine-tuned to cancel a large portion of the S-D contribution generated by the tensor part of the OPE. Further, the inclusion of the O(p 2 ) counter terms D d and D f between the S-wave elastic channels in addition to the D SD term results only in a minor change in the fit (see the red solid curves for Scheme F) with the χ 2 /d.o.f = 1.10. Finally, by adding also the OEE interaction (Scheme G) one obtains results which lie on top of those for Scheme F in Fig. 4 and which are therefore not shown.
We are now in a position to re-analyse the net effect from the OPE within Scheme F, which allows one to achieve an overall good description of the currently available experimental data with χ 2 /d.o.f around unity. In Fig. 5, we compare the results of the fits for Scheme F (red solid curves) with those where pions are switched off (Scheme H -the black dashed curves) to conclude that the intermediate momentum scale p typ generated by the pionic dynamics can be absorbed to a very large extent into the re-definition of the O(p 2 ) contact terms (essentially D SD ). Thus, the current data in the 1 +− channel basically do not leave room for the residual long-range chiral dynamics from the OPE. Meanwhile, the difference between the results for Schemes H (black dashed curves) and A (blue dotted curves), which is best seen in the BB * channel, represents the net effect of the O(p 2 ) contact interactions on the line shapes. Finally, in Fig. 6 we present the results for Scheme G including the uncertainties which correspond to a 1σ deviation in the parameters including correlations.
In summary, the results presented in this section demonstrate that (i) the data can be equally well described with or without the central S-wave OPE potential since the latter can be absorbed into a re-definition of the S-wave shortrange contact terms;  (ii) the inclusion of D waves affects the line shapes noticeably which confirms the claims found in the literature that D waves are important in the near-threshold charmoniumlike and bottomonium-like systems [44,61]. However, the current data call for the promotion of the O(p 2 ) S-wave-to-D-wave counter term D SD to lower order and for tuning this term to balance the S-D dynamics from the OPE. Thus, given the accuracy of the present experimental data (especially those in the BB * channel), the residual effect from the OPE on the line shapes is small; (iii) the effect of the OEE interaction is negligibly small; (iv) the data are essentially consistent with HQSS constraints imposed on the potential; (v) no indication for the importance of O(p 2 ) contact interactions in the inelastic channels is seen in the data. (a)

Dependence of the results on the regulator
In this chapter, we investigate if the promotion of the S-D counter term is called for by the renormalisation of the leading order amplitudes. In Fig. 7 we show the regulator dependence of the results corresponding to the fits F and C, that is, those fits with and without the O(p 2 ) contact terms in addition to the OPE. The results for fit C show a clearly visible cut-off dependence in the BB * channel when the (sharp) cut-off in the LSEs, treated as a hard scale, is varied in a reasonable range from 800 to 1300 MeV (cf. black long-dashed, blue dashed and green dotted-dashed curves). Note that in an effective field theory where a potential, which is then iterated in a scattering equation, is expanded in terms of a given expansion parameter one cannot expect a complete regulator independence of observable quantities at a given order. This implies that the regulator should not be chosen too large and that regulator effects of sub-leading order are common. It is thus difficult to judge, if the mentioned variations call for an additional counter term at LO or not. However, it is certainly interesting to observe that for the same cut-off variation the curves in fit F remain unchanged indicating that the theory is fully renormalised. The same pattern is observed in the other channels which are therefore not shown. The red curves in Fig. 7 (both for the fits F and C) correspond to a much harder cut-off, Λ = 550 MeV, which leaves basically no room for the contributions generated by the coupled-channel dynamics from the OPE governed by the scale p typ ∼ 500 MeV. It is, therefore, not surprising that the red curves for the two fits are almost identical. It is also instructive to note that for the fit F all theoretical curves, including the one for the hardest cut-off, provide an almost equally good description of the data. This should have been expected indeed since the D SD term basically absorbs all the dynamics related to the scale p typ . Therefore, the spread with the cut-off in the results for the fit F is minor in the whole interval of considered cut-offs. 5 The pole positions and the nature of the Z b states In this section we discuss the extraction of the poles of the amplitude in the complex plane. In general, in order to search for the poles in a multi-channel problem a multi-sheet Riemann surface in the complex energy plane needs to be invoked. However, we consider the four-sheet Riemann surface corresponding to the two elastic channels only, because all inelastic thresholds are remote and their impact on the poles of interest, which are located near the elastic thresholds, is expected to be minor. 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 variable which is traditionally denoted as ω [75] and which, for a given energy E = M − m B − m B * , is defined via 1) where µ 1 and µ 2 are the reduced masses in the first (BB * ) and in the second (B * B * ) elastic channel, respectively. 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 reads 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.
The corresponding regions in the ω-plane are depicted in the first plot in Fig. 8. The thick solid line corresponds to the real values of the energy lying on RS-I. It is easy to see that the physical region between the two thresholds corresponds to |ω| = 1, with both Re(ω) and Im(ω) positive, and the thresholds at E = 0 and E = δ are mapped to the points ω = ±i and ω = ±1, respectively. The inelastic channels in the energy plane are interpreted as one additional effective remote channel with the momentum k in and, in order to find all relevant poles, both possibilities with Im k in > 0 and Im k in < 0 are considered. 6 As it is now effectively a three-channel problem, the number of Riemann sheets doubles and so does the number of poles representing the physics for each state. Hence, each pole has now its mirror partner located on a different sheet of the eight-sheeted Riemann surface but corresponding to the same physical state. The pole positions in the ω-plane for the fit Schemes A and G from the previous section are shown in the right panel of Fig. 8. As long as the inelastic channels are switched off, one is back to the two-channel problem with one relevant pole corresponding to each state. For example, the poles for the fit Scheme A but without inelastic channels are shown by the crosses (x) in Fig. 8: the cross on the imaginary axis on RS-II (close to ω = i) corresponds to a virtual BB * state pole associated with the Z b state while the other pole, residing on RS-IV (close to ω = 1), represents the physics related to the Z b -it is a virtual state in the B * B * channel slightly shifted to the complex plane due to coupled-channel effects. 7 When the effective inelastic channel is on, one arrives at a pair of poles for each state: one for Im k in > 0 and one for Im k in < 0 -the corresponding solutions are labeled by "+" and "−", respectively. For the Z b , the two resulting poles on RS-II + and RS-II − are symmetric with respect to the imaginary axis in the omega plane which results in complex-conjugate solutions for the energies (see Table 4). Unlike the Z b , the poles for the Z b on RS-III − and RS-IV + are slightly asymmetric due to coupled-channel effects. In any case, it is apparent that the poles obtained in the full calculation, including inelastic channels, reside in the vicinity of those found without inelastic channels, that implies that the role played by the inelastic channels is sub-leading in line with a molecular interpretation of the Z b states.
The corresponding energies evaluated relative to the relevant elastic threshold, 2) 6 In the actual calculation, the various inelastic channels have different momenta. However, when searching for the poles, all inelastic channels are assumed to be on the same sheet, that is the imaginary parts of all inelastic momenta are synchronised to be either positive or negative. 7 Both these poles have counterparts on RS-IV around ω = −i and ω = −1, respectively, that is far away from the physical region. Analogous additional poles also appear for the other cases discussed below, but will not be mentioned anymore because they do not affect the line shapes.  are listed in Table 4 and are visualised in Figs. 8 and 9 (the errors in the pole positions correspond to a 1σ deviation for the whole parameter list). Note that the sign convention is such that positive energies refer to above-thresholds poles -see definition (5.2).  Table 4. The energies E Z b and E Z b (see the definition in Eq. (5.2)) for the fit Schemes A and G. The energies denoted as A − and G − stand for the poles for the fit Schemes A and G, respectively, with all the inelastic channels on their unphysical (Im k in < 0) Riemann sheets. The energies denoted as A + and G + are for the poles for the fit Schemes A and G, respectively, with all the inelastic channels on their physical (Im k in > 0) Riemann sheets. The errors correspond to a 1σ deviation in all fitted parameters.
The two poles for Fit A representing the results based on S-wave contact interactions (shown as an up-(red) and down-pointing (green) triangles in Figs. 8 and 9) are consistent within 1σ with those obtained in Ref. [27] -see Fig. 8 of the quoted paper. As one can see from Table 4 Tables 2 and 3), the overall description of the data in our best fit G is only slightly better than in fit A, as one can see from the line shapes shown in Fig. 5 and from the similar values of χ 2 /d.o.f. quoted in Table 2. Thus, both fits describe the Z b state as a shallow virtual state located below the BB * threshold. Meanwhile, the Z b state is consistent with both a virtual state and an above-threshold resonance interpretation, with the latter option preferred by the fit G.

Summary
This work continues a series of papers aimed at a systematic description of the line shapes of near-threshold resonances in general and, in particular, at understanding the nature and properties of the Z ± b (10610) and Z ± b (10650) states in the spectrum of bottomonium. Unlike the previous papers in the series [26,27], in this work, we do not resort to an analytic parameterisation for the line shapes but rely on an EFT approach to calculate the line shapes explicitly. In particular, in order to obtain the production and scattering amplitudes, we construct the effective potential and iterate it to all orders using the coupledchannel Lippmann-Schwinger equations numerically. This allows us to verify the accuracy of the practical parameterisation suggested in Ref. [26] and to estimate the role played by the π-and η-exchanges, which cannot be straightforwardly incorporated into the scheme of Refs. [26,27] because of their non-separable form.
The results of this work can be summarised as follows: • We find that the distributions obtained from the direct numerical solution of the LSEs with just S-wave contact potentials provide a nearly identical description of the data to that achieved using the parameterisation derived in Refs. [26,27,66].
• We include the OPE interaction on top of the contact potentials in the elastic channels and demonstrate that the inclusion of the S-wave OPE affects the line shapes only marginally. After a re-fit to the data required to appropriately renormalise the shortrange interactions, basically the whole S-wave OPE contribution can be absorbed to a re-definition of the S-wave contact interactions leaving only a small residual effect. Therefore, we do not support the claim of Ref. [64] that OPE changes the line shape of near threshold states by about 30%.
• For the actual value of the splitting of the two elastic thresholds, δ ≈ 45 MeV, the momentum scale p typ 500 MeV which controls the coupled-channel (BB * -B * B * ) dynamics is relatively large -see Eq. (1.1). For such momenta, the role played by the OPE in D waves turns out to be significant even after a re-fit, so that one cannot neglect this effect in order to extract the resonance parameters with a sufficiently high accuracy. However, the resulting significant distortion of the line shapes from the OPE is not supported by the currently available experimental data in the BB * channel: Indeed, a clear bump structure around the B * B * threshold unavoidably generated by the S-D OPE transitions is not seen in the most recent data. In order to cure this, we include additional contact terms allowed by chiral symmetry at order O(p 2 ) to find that the resulting line shapes are in a very good agreement with the data, with a reduced χ 2 /d.o.f. close to 1. The role played by various O(p 2 ) contact interactions is different: on the one hand, a single contact term which contributes to the S-to-D-wave transitions absorbs the largest part of the S-D OPE piece and, in addition, brings an additional residual contribution -both effects together improve the quality of the fit considerably; on the other hand, two S-S contact interactions play a sub-leading role resulting only in a marginal change in the fits. Finally, we conclude that after the inclusion of the O(p 2 ) contact interactions the residual effect from the OPE is marginal.
• There are indications that the inclusion of the S-D contact term at leading order is required by renormalisation. However, to make this statement more sound a complete calculation to next-to-leading order would be necessary. Since this would come with additional parameters that need to be fixed by data, we postpone this effort until improved data become available.
• We find that the data are consistent with HQSS symmetry constraints imposed on the potential.
• The data do not call for the inclusion of O(p 2 ) contact interactions in the elastic-toinelastic transitions.
• We extract the position of the poles responsible for the Z ± b (10610) and Z ± b (10650) and find them to reside on the unphysical Riemann sheets just below (Z b ) or just above (Z b ) the corresponding elastic threshold.
Before closing, we would like to stress that the observed almost complete cancellation between the OPE and the additional O(p 2 ) S-D contact terms is very puzzling. The pion plays a very special role in QCD due to its intimate connection to the spontaneous breaking of chiral symmetry. In addition, since the D * Dπ coupling is known (see Eq. (2.16)), the OPE comes without adjustable parameters. Still, data demand an almost complete cancellation of this part of the potential. Whether or not this cancellation is accidental in the system at hand or if it has deeper reasons in QCD calls for further studies which, however, go beyond the scope of this work. Clearly, more accurate experimental data for the Z b 's, especially in the elastic channels, and, hopefully, new data for their spin partner states would be of great relevance and importance for such studies. Sino-German CRC 110 "Symmetries and the Emergence of Structure in QCD". Work of V.B. and A.N. was supported by the Russian Science Foundation (Grant No. 18-12-00226).

A Effective Lagrangians
The effective Lagrangian describing isovector B ( * )B( * ) scattering at low energies reads where the terms in the first row in Eq. (A.1) stand for the leading heavy hadron chiral perturbation theory Lagrangian of Refs. [77][78][79], written in the two-component notation of Ref. [80]. The terms in the second row represent the Lagrangian for anti-heavy mesons. The heavy mesons and anti-heavy mesons interact via the remaining terms in the Lagrangian: the third row in Eq. (A.1) corresponds to O(p 0 ) S-wave contact interactions [38,67], while the last three rows represent the O(p 2 ) contact terms with two derivatives. The contact terms ∝ D 10 and D 11 contribute to S-wave interactions while the last term ∝ D 12 is projected out to give rise to the S-D transitions. Since we are only interested in the S-S and S-D transitions for the B ( * )B( * ) scattering, all the terms of the kind ∝ ∇ i H † ∇ j H contributing to P waves were dropped. In Eq. (A.1) H a = P a + V i a σ i stands for the heavy meson superfield combining a pseudoscalar (P a ) and a vector (V a ) meson, a and b are SU (2) isospin indices, and the isospin matrices are normalised via the trace as τ A ab τ B ba = 2δ AB . The charge and Hermitian conjugate operators for the superfield H a read and they transform as  where x = cos(pq) and tiny corrections ∝ δ p 2 typ /2m B were neglected. It is the term ∝ x that eventually supports higher partial waves in the elastic channel. Then the D waves are suppressed relative to the S waves by the factor (B. 3) The values of q range from about 200 MeV, for the transition to the Υ(3S)π channel, to 1.1 GeV, for the transition to the Υ(1S)π channel. It is easy to verify that, for such values of q, the second factor in Eq. (B.3) is close to unity, that provides a justification for the estimate used in the main text -see the discussion above Eq. (2.4).

C Flavour symmetry and the pseudoscalar exchange potentials
In this appendix we discuss the pseudoscalar exchange potential between the heavy mesons. For definiteness, we stick to the BB * channel as to an example. The flavour projector for a pseudoscalar (P orP ) and a vector (V or V ) for the 1 +− quantum numbers reads where the overall factor 1/2 ensures the proper normalisation of the projectors, Here, like in Appendix A, the standard normalisation for the isospin matrices was used, Tr[τ A τ B ] = 2δ AB .
Then, the flavour factor involving isospin and C-parity for the OPE diagram can be evaluated as 1 4 Tr V † τ A P † −P † τ A V † P † τ a V V τ BP − P τ BV V † τ aP where at the last step above an easily verified relation Tr τ A τ a τ B τ a = −2δ AB (C.4) was used. It is this part of the calculation which makes difference between the OPE and OEE potentials. Indeed, in the latter case the trace takes the form Tr τ A1 τ B1 = 2δ AB , (C.5) where the unit matrices correspond to the η emission/absorption vertices which substitute the Pauli matrices from the pion vertices -see Eq. (C.4). Finally, considering an additional factor 1/ √ 3 which comes from the 8-th Gell-Mann matrix for the SU (3) group -see Eq. (2.13) -and which, therefore, enters each η-vertex, one arrives at the following list of changes needed to proceed from the OPE potential to the OEE one: (i) the flavour factor +1, for the OPE in the isovector channel, should be replaced by −1, for the OEE; (ii) the pion mass m π should be replaced by the η mass m η , and (iii) the pion coupling constant g b should be replaced by the η coupling constant g b / √ 3. Interestingly, if the flavour symmetry group SU (3) is extended to U (3) then the potentials from the π-, η-, and η -exchange cancel each other identically in the strict flavour symmetry limit [65]. However, because of the U (1) anomaly, there are no reasons to expect the η to possess properties of the Goldstone boson [81][82][83][84].