The $Z_b$ structures in a constituent quark model coupled-channels calculation

The $Z_b(10610)^\pm$ and $Z_b(10650)^\pm$ are two bottomonium-like structures discovered in the $\pi h_b(mP)$, $\pi \Upsilon(nS)$ and $B^\ast\bar B^{(\ast)}+h.c.$ invariant mass spectra, where $m=\{1,2\}$ and $n=\{1,2,3\}$. Their nature is puzzling due to their charge, which forces its minimal quark content to be $b\bar b u\bar d$ ($b\bar b d\bar u$). Thus, it is necessary to explore four-quark systems in order to understand their inner structure. Additionally, their strong coupling to channels such as $\pi \Upsilon$ and the closeness of their mass to $B^\ast\bar B^{(\ast)}$-thresholds stimulates a molecular interpretation. Within the framework of a constituent quark model which satisfactorily describes a wide range of properties of (non-)conventional hadrons containing heavy quarks, we perform a coupled-channels calculation of the $I^G(J^{PC})=1^+(1^{+-})$ hidden-bottom sector including $B^{(\ast)}\bar B^{\ast}+h.c.$, $\pi h_b$, $\pi \Upsilon$ and $\rho\eta_b$ channels. We analyze the line shapes in the different channels, describing the $\Upsilon(5S)\to \pi B^{(*)}\bar B^{(*)}$ by means of the $^3P_0$ model. Since our description of the line shapes promising, we perform the same coupled-channels calculation for the $Z_b$'s with $J^{--}$, where $J=\{0,1,2\}$. This allows us to obtain a fair description of the corresponding line shapes. The study of the analytic structure of the $S$-matrix suggests that the experimental $Z_b$ structures arise as a combination of several poles with $J^{PC}=0^{--}$, $1^{\pm-}$ and $2^{--}$ quantum numbers nearby the $B\bar B^\ast$ and $B^\ast\bar B^\ast$ thresholds.


I. INTRODUCTION
During almost two decades, tens of charmonium-and bottomonium-like states, the so-called XYZ mesons, have been observed at B-factories (BaBar, Belle and CLEO), τ -charm facilities (CLEO-c and BESIII) and also proton-(anti)proton colliders (CDF, D0, LHCb, ATLAS and CMS). The common characteristic of most of these states is that, although their quantum numbers are compatible with naive QQ states (Q either c-or b-quark), their masses, static properties and decay patterns point out to more complex structures involving higher Fock-state components.
It was not until 2011 when undeniable evidences of exotic mesons with forbidden quantum numbers for a quark-antiquark pair were observed. The charged bottomonium-like states Z b (10610) ± and Z b (10650) ± were identified by the Belle Collaboration [1] as peaks in the invariant mass distribution of the π ± h b (mP ) (m = 1, 2) and π ± Υ(nS) (n = 1, 2, 3) subsystems when the Υ(10860) resonance decays into two pions plus an h b or Υ. Shortly after, the Belle Collaboration [2,3] confirmed their existence when investigating the elastic BB * and B * B * channels. Moreover, the quantum numbers of the Z b 's were analyzed to be I G (J P C ) = 1 + (1 +− ) [4] and so it appears in the Review of Particle Physics (RPP) of the Particle Data Group (PDG) [5]. Note herein that a 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 flavor symmetry (HQFS) should be very accurate for these systems. In fact, one can predict their HQFS partners in the cc sector and assign them to the Z c (3900) ± /Z c (3885) ± and Z c (4020) ± charmonium-like structures discovered in the πJ/ψ, πh c and D * D( * ) + h.c. invariant mass spectra [6,7]. Their strong coupling to channels such as πJ/ψ and the closeness of their mass to D * D( * ) -thresholds stimulates a molecular interpretation and, among an extensive literature concentrated in these charmonium-like states, we performed in Ref. [8] a coupled-channels calculation of the I G (J P C ) = 1 + (1 +− ) sector including D ( * )D * + h.c., πJ/ψ and ρη c channels in the framework of a constituent quark model. Therein, we observed that the mesonmeson interactions are dominated by the non-diagonal πJ/ψ − D * D( * ) and ρη c − D * D( * ) couplings indicating that the Z c (3900) ± /Z c (3885) ± and Z c (4020) ± are unusual states. The study of the analytic structure of the S-matrix allowed us to conclude that the point-wise behavior of the line shapes in the πJ/ψ and DD * invariant mass distributions is due to the presence of two virtual states that produce the Z c peaks. This calculation was recently extended to the open-strange sector with the description of the Z cs (3985) − [9], with ccsū minimal quark content.
Turning our attention back to the Z b states, both tetraquark [10][11][12] and molecular [13][14][15][16][17][18][19][20][21][22] pictures are claimed to be in fair agreement with the experimental arXiv:2107.02544v1 [hep-ph] 6 Jul 2021 data. However, the Z b (10610) ± and Z b (10650) ± exotic states are located quite close to the BB * and B * B * thresholds, respectively; and, moreover, these channels are the dominant decay channels of the Z b 's. This provides a strong support for their molecular interpretation. Reference [23] argued that the Z b (10610) ± and Z b (10650) ± could be simple kinematical cusps; however, general arguments [24] and explicit calculations [25,26] demonstrate that narrow bumps as the Z b 's signals in the BB * and B * B * channels necessitate near-threshold poles.
Lattice QCD studies of the Z b 's are scarce because the I G (J P C ) = 1 + (1 +− ) hidden-bottom sector presents severe challenges. If the Lüscher method was used to study the poles of the Scattering matrix one should include, at least, seven coupled two-meson channels. Furthermore, the original Lüscher approach for two-particle scattering is not valid above the three-particle threshold; though its extension to the three-body scattering processes is underway [50,51]. Despite of these difficulties, preliminary lattice studies have been reported in Refs. [52,53]. Moreover, a Born-Oppenheimer approximation [54,55], inspired by the study of these systems in [52,53], has been recently applied to the Z b states within the framework of lattice QCD [56,57].
The structure of the present manuscript is organized as follow. After this introduction, the theoretical framework is briefly discussed in Sec. II. Our analysis of the obtained results: masses, decay widths and line shapes, can be found in Section III. Finally, we summarize and give some conclusions in Sec. IV.

II. THEORETICAL FORMALISM
A. Constituent quark model In the absence of current quark masses, Quantum Chromodynamics (QCD) is invariant under chiral transformations. However, this symmetry is not realized in Nature indicating that it is spontaneously broken in QCD. Among its many consequences, the most outstanding are the development of a constituent quark mass which depends on the quark momentum, M = M (q 2 ) with M (q 2 = ∞) = m q being the current quark mass that appears in the QCD Lagrangian, and the existence of Goldstone-bosons whose exchanges between the light quarks produce interactions among them.
We use a so-called constituent quark model (CQM) which mimics the mentioned phenomena based on the following effective Lagrangian at low-energy [78]: where U γ5 = e iλaφ a γ5/fπ is the matrix of Goldstoneboson fields that can be expanded as follows in which the first term gives rise to the constituent quark mass, the second one describes the pseudoscalar-meson exchange interaction among quarks and the main contribution of the third term can be modeled by means of a scalar-meson exchange potential. Another well-known non-perturbative phenomenon is the confinement of quarks inside hadrons. Lattice-QCD studies have demonstrated that multi-gluon exchanges produce an attractive linearly rising potential proportional to the distance between infinite-heavy quarks [79]. However, the spontaneous creation of light-quark pairs from the QCD vacuum may give rise at the same scale to a breakup of the color flux-tube [79]. We have tried to mimic these two phenomenological observations by the following expression: Here, a c and µ c are model parameters, and the SU(3) color Gell-Mann matrices are denoted as λ c . Finally, one expects that the dynamics of the boundstate system is governed by QCD perturbative effects at short inter-quark distances. This is taken into account by the one-gluon exchange potential derived from the following vertex Lagrangian with α s an effective scale-dependent strong coupling constant which allows us to consistently describe light, strange and heavy quark sectors, and whose explicit expression can be found in, e.g., Ref. [80]. A detailed physical background of the quark model can be found in Refs. [58,59]. The model parameters and explicit expressions for the potentials can be also found therein. We want to highlight here that the interaction terms between light-light, light-heavy and heavy-heavy quarks are not the same in our formalism, i.e. while Goldstone-boson exchanges are considered when the two quarks are light, they do not appear in the other two configurations: light-heavy and heavy-heavy; however, the one-gluon exchange and confining potentials are flavorblind.

B. Resonating group method
The Resonating Group Method (RGM) [81] allows us to describe the strong interaction at the meson level. We are going to consider mesons as quark-antiquark clusters and the effective cluster-cluster interaction emerges from the aforementioned microscopic interaction among constituent quarks.
The wave function of a system composed of two mesons, A and B, with distinguishable quarks, can be written as 1 where, e.g., φ A ( p A ) is the wave function of the meson A with p A the relative momentum between its quark and antiquark. The relative motion of the two mesons is taken into account by the wave function χ α ( P ). The projected Schrödinger equation for the relative wave function can be written as follows: where E is the energy of the meson-meson system. The direct potential, RGM V αα D ( P , P i ), can be written as whereas RGM V αα R ( P , P i ) is the quark re-arrangement potential and represents a natural way to connect mesonmeson channels with different quark content such as πΥ and BB * . It is given by where P mn is an operator that exchanges quarks between clusters.
The Rayleigh-Ritz variational principle is used to determine the eigen-values and -functions of the two-body Schrödinger equation that describes mesons. This is because its simplicity and flexibility; however, the choice of basis to expand the wave function solution is of great importance. We use the Gaussian Expansion Method [82] which provides enough accuracy and simplifies the subsequent evaluation of the needed matrix elements. With the aim of optimizing the Gaussian ranges employing a reduced number of free parameters, Gaussian trial functions, whose ranges are given by a geometrical progression, are introduced [82]. This choice produces a dense distribution at short distances enabling better description of the dynamics mediated by short range potentials.
The coupled-channels RGM equation, Eq. (6), can be rewritten as a set of coupled Lippmann-Schwinger equations of the following form where α labels the set of quantum numbers needed to uniquely define a certain partial wave, V α α (p , p) is the projected potential that contains the direct and rearrangement potentials, and E α (p ) is the energy corresponding to a momentum p , written in the nonrelativistic case as: Herein, µ α is the reduced mass of the meson-meson system corresponding to the channel α, and ∆M α is the difference between the meson-meson threshold and the one we shall take as a reference. The coupled-channels Lippmann-Schwinger equation is solved using a generalization of the matrix-inversion method [83] in order to include channels with different thresholds. Once the T -matrix is calculated, we determine the on-shell part which is directly related to the scattering matrix (in the case of non-relativistic kinematics): with k α the on-shell momentum for channel α.
Since our aim is to explore within the same formalism the existence of states below and above meson-meson thresholds. All the potentials and kernels shall be analytically continue for complex momenta; this allows to find the poles of the T -matrix in any possible Riemann sheet.

C. Three-body production amplitude
The experimental line shapes analyzed in this manuscript are calculated in e + e − collisions at the center of mass energy of the Υ(10860) meson. Thus, it is reasonable to think that a large portion of the Z b candidates emerges from mechanisms in which the Υ(10860) dessintegrates into a pion plus a Z b .
In order to shed some light into the Z b decay line shapes and to analyze possible contribution from alternative J P C sectors, it is relevant to properly account for the momentum dependence of the Z b s production vertex, as it would be questionable to assume a constant production vertex for channels in a relative P wave, as the πh b (mP ) channels in the J P C = 1 +− sector.
In a molecular interpretation for the Z b candidates, the reaction Υ → π(M 1 M 2 ) involves the creation of two quark-antiquark pairs (see Fig. 1), where M 1 M 2 is the final two-body state of the coupled-channels calculation, i.e., M 1 M 2 = {πΥ, πh b , BB * , B * B * }. A simple approach for such reaction can be obtained from the quark-antiquark pair creation 3 P 0 model [84][85][86]. The transition operator for the aforementioned reaction can be written as, where µ i (ν i =μ i ) are the quark (antiquark) quantum numbers for the pair i = {1, 2}, and γ = 2 5/2 √ 2π 1/2 γ with γ = g 2m is a dimensionless constant that gives the strength of the qq pair creation from the vacuum. The subscript brakets in the operator represent the {C, I, S, J} quantum numbers of the created qq pair.
From this transition operator we can obtain the potential V3 P0 (k), where k is the relative momentum of the final two mesons in the coupled-channels calculation. For this potential, we employ the same meson wave functions obtained from the CQM model, as done for the coupled-channels calculation. A relevant key of this approach for the Z b production is that the 3 P 0 mechanism only allows the Υ(10860) → π(B ( * )B( * ) ) transitions, so the production of hidden-bottom channels such as πΥ(nS) or πh b (mP ) only arise through a further quark rearrangement process, i.e., B ( * )B( * ) → πΥ(nS) or B ( * )B( * ) → πh b (mP ).

D. Production line shapes
The Z b → M 1 M 2 production is described diagramatically as in Fig. 2, and it can be written in terms of the following Lorentz-invariant production amplitude: where β ( ) denotes a given M 1 M 2 channel in the coupledchannels calculation and the amplitude is summed up over the different J P C considered in this work. In order to explore the description of line shapes, we add the amplitudes A β and phases θ β to take into account different production weights for each channel M 1 M 2 . It is worth reminding that the background contribution only affects the BB * and B * B * channels, as the hidden-bottom channels are OZI-suppressed. A resonance in its rest frame, which decays into a three-body final state π(M 1 M 2 ), has a differential decay rate given by where M is the Lorentz-invariant production amplitude from the 3 P 0 amplitude described above, defined in Eq. (13), and dΦ is the three-body phase-space which can be written as Therefore, we can express dΓ in terms of the invariant mass spectrum of the two-meson channel, m M1M2 , as where (k M1M2 , Ω M1M2 ) is the on-shell momentum of the M 1 M 2 pair and Ω πZ b is the solid angle of the π in the center-of-mass reference frame of π(M 1 M 2 ) at energy √ s. The on-shell momenta are given by , with m ± = m 1 ± m 2 . Integrating the angles, we have with β the quantum numbers of the channel M 1 M 2 and M β is averaged over Υ(5S) spin states.
Experimentalists actually measure events in the σ(e + e − → π ± Z ∓ b ) × B (Z ∓ b → M 1 M 2 ) process; to describe the data, we need to add a normalization factor to translate the decay rate into events: In order to fit {A M1M2 , θ M1M2 , N M1M2 } for each channel, we minimize a global χ 2 function using all the available experimental data for channels πΥ, πh b , BB * and B * B * : The theoretical lineshapes are, then, convolved with a σ = 6 MeV Gaussian, to take into account the experimental resolution [1][2][3].

III. RESULTS
It is worth emphasizing firstly that Heavy Flavour Symmetry (HFS) dictates the B ( * )B * and D ( * )D * interactions are practically identical. This happens because the B ( * ) and D ( * ) wave functions are similar and, regarding the meson-meson interaction, the heavy quark c, or  b, acts as an spectator and the dynamics is govern by the light quark pair. However, the available kinetic energy of the involved channels is reduced in the bottom sector due to the larger mass of the b-quark, favoring the creation of bound states. We begin performing a coupled-channels calculation for the I G (J P C ) = 1 + (1 +− ) sector, including the fol- In fact, the binding-energy contribution of this channel would come from the 1/m q mq suppressed spin-orbit term, which would be able to flip the spin of the heavyquark pair.
The first column of Tables I and II show our findings for the I G (J P C ) = 1 + (1 +− ) hidden-bottom sector. Contrary to the Z c case [8], the B ( * )B * interaction is strong enough to produce two real bound states below the BB * and the B * B * thresholds. Tables I and II collect the largest meson-meson component in the wave function of our candidates for the Z b (10610) ± and Z b (10650) ± states; one can see that the dominant wave-function's component is BB * and B * B * , respectively, which facilitates their interpretation as B ( * )B * molecular states.
Our prediction for the Z b 's total decay widths, broken down into the different meson-meson final state contributions, are also presented in Tables I and II. We are not able to reproduce the total decay widths reported by the Belle collaboration for neither Z b (10610) ± nor Z b (10650) ± . Moreover, their coupling to the πh b channel is very little (not collected in the Table) with the Z b 's interpreted as molecular states of B * B * with quantum numbers I G (J P C ) = 1 + (1 +− ).
With this theoretical description of the Z b 's at hand, we computed their production rates and compare them  with the experimental data. We consider that the Z b (10610) ± and Z b (10650) ± states are coming from a Υ(5S) → π(B ( * )B( * ) ) → π(M 1 M 2 ) vertex described by means of a 3 P 0 amplitude with a center-of-mass energy √ s = 10.865 GeV, according to Refs. [1][2][3]87]. On a first step, we can evaluate the goodness of the proposed 3 P 0 production amplitude by analyzing the decay branching ratios of Υ(5S) → πB ( * )B( * ) , integrating Eq. (18). Experimental values for such branchings have been measured in Ref. [88]: over a total width of Γ exp Υ = 37 ± 4 MeV [5]. The only parameter of the production amplitude is the value of γ of the 3 P 0 model. We will take the same value fixed for two-body decays of bottomonium studied in Ref. [59], considering a 10% uncertainty, viz. γ = 0.205 ± 0.020. Without considering rescattering contributions in Eq. (18), such value of γ gives, which are in good agreement with the experimental values. If we analyze the contribution of each B ( * )B( * ) partial wave to the total (Table III), we see that the S and P waves are of the same order, so we expect a similar contribution of S and P waves in the line shapes. Hence, as the production through a P wave is significant, in order to deepen our understanding of the Z b 's  inner structure, we have performed three further coupledchannels computations for the quantum numbers J −− , with J = {0, 1, 2}. Indeed, this allows us to include the πh b (mP ) (m = 1, 2) channels in our calculation in a relative S wave in the 1 −− sector 2 . For completeness, we add the BB channel in the coupled-channels calculation for those sectors.
One can see in Tables I and II the S-matrix pole structure for the J −− hidden-bottom sector, also summarized in Fig. 3. As the B ( * ) B * can only be in a relative P -wave for the J −− , we mostly obtain virtual states below the B * B( * ) thresholds. A near-threshold resonance above B * B * channel is found for the 1 −− sector, with mass close to the Z b (10650) ± state. All the states predicted have a dominant B ( * ) B * component in their wave functions and thus their interpretation as molecular states is presumably.
The computed decay widths of the three states, and their decomposition into the different meson-meson channels, are also collected in Tables I and II. We are predicting the same order of magnitude, slightly higher, for the total decay widths reported experimentally: (18.4 ± 2.4) MeV for the Z b (10610) ± , and (11.5 ± 2.2) MeV for the Z b (10650) ± [5]. Moreover, the most important decay channels are πΥ(nS) with n = 1, 2, 3 for the Z b (10610) ± assignments, and B * B * for the three candidates of the Z b (10650) ± . Let us finally mention here that the distinct states around the B * B * and BB * thresholds with different quantum numbers J P C could be being interpreted as single Z b (10610) ± and Z b (10650) ± resonances in the experiments.
The point-wise behavior of the line shapes for an invariant mass between 10.55 and 10.75 GeV can be seen in Fig. 4 for the BB * and B * B * final decay channels, and in Fig. 5 for the πΥ(nS) (n = 1, 2, 3) and πh b (mP ) (m = 1, 2) final decay channels. We find a good agreement with the B ( * )B * experimental data in Fig. 4, due to the production from all sectors considered in this work. The J −− poles close to the experimental Z b (10610) ± and Z b (10650) ± thresholds influence the BB * and B * B * line shapes, respectively. The Z b (10650) ± peak in the BB * line shape is seen as a small bump around 10.65 GeV, smeared by the detector resolution.
The description of the line shapes for hidden-bottom channels πΥ(nS) and πh b (mP ) is shown in Fig. 5. It is worth remarking that all these channels are produced solely through rearrangement potentials from the B ( * )B( * ) channels. The πΥ(nS) (n = 1, 2, 3) channels show two bumps which can be assigned to the two Z b 's structures. A general good agreement is found for these line shapes, though the intensity of the peaks for the πΥ(3S) channel is somehow lower than the experimental data.
Regarding the πh b (1P ) and πh b (2P ) line shapes, Figure 5, the agreement to the experimental data is fair, showing a dominance of the second peak, which mostly originates from the 1 −− sector. The line shape shows a nearly constant production with the invariant mass until the B * B * mass threshold and then it falls-off quickly to zero, describing the experimental Z b (10650) ± structure. The first peak, corresponding to the Z b (10610) ± structure is, however, smeared, which points to a weaker BB * → πh b (mP ) interaction in the 1 +− sector compared to the 1 −− B * B * → πh b (mP ).  [3], whose statistical uncertainty is represented by a vertical line.

IV. SUMMARY
The Z c 's are very peculiar states, different from other potential molecular candidates of the hidden-charm spectrum. In order to clarify their inner structure, within the framework of a constituent quark model that satisfactorily describes a wide range of properties of (non-)conventional hadrons containing heavy quarks, we performed in Ref. [8] a coupled-channels calculation of the I G (J P C ) = 1 + (1 +− ) ccud (ccdū) sector around the mass energies of the Z c (3900) ± /Z c (3885) ± and Z c (4020) signals. We found that the diagonal interaction between the D ( * )D( * ) channels is, too, suppressed to develop resonances, being the non-diagonal rearrangement interaction due to the coupling with other meson-meson channels responsible for their virtual-state nature, needed to describe the relevant line shapes. The same approach has FIG. 5. Line shapes of the πΥ(nS) (n = 1, 2, 3) and πh b (mP ) (m = 1, 2) channels (see legend in the abscissa axis). Same legend as Fig. 4. The black solid points are the experimental data from Refs. [1,2,87], whose statistical uncertainty is represented by a vertical line.
also been applied to the ccsū sector pursuing a description of the Z cs (3985) − [9].
We have naturally extended herein our analysis to the the Z b (10610) ± and Z b (10650) ± signals observed in 2011 by the Belle Collaboration, viz. a coupled-channels calculation of the I G (J P C ) = 1 + (1 +− ) bbud (bbdū) sector, including B ( * )B * , πΥ(nS) (n = 1, 2, 3), πh b (mP ) (m = 1, 2) and ρη b channels, has been performed within the same constituent quark model framework. In this initial calculation, we found that the πh b channel coupling to BB * was rather small (due to the supressed spin-orbit term), and zero for the B * B * channel when only the quantum numbers I G (J P C ) = 1 + (1 +− ) are available for the Z b (10610) ± and Z b (10650) ± states. Thus, we per-form the same coupled-channels calculation but considering too the following quantum numbers: J P C = 0 −− , 1 −− and 2 −− for the Z b 's. This preference is based on two main reasons: On the one hand, it allows us to include the πh b final states in S-wave, with a non-zero coupling to the B * B * channel. On the other hand, the description of the Υ(10860) → π ± Z ∓ b vertex using the 3 P 0 model points to a significant contribution from P -wave B ( * )B( * ) channels.
Several S-matrix poles were identified around the BB * and B * B * thresholds, with masses and widths compatible with the Z b (10610) ± and Z b (10650) ± resonances, and we obtained a fair description of the corresponding line shapes due to the combination of the production amplitudes of all those quantum numbers.