The charged $Z_c$ and $Z_b$ structures in a constituent quark model approach

The nature of the recently discovered $Z_c$ and $Z_b$ structures is intriguing. Their charge forces its minimal quark content to be $Q\bar Q q\bar q$ (where $Q=\{c,b\}$ and $q=\{u,d\}$). In this work we perform a molecular coupled-channels calculation of the $I^G(J^{PC})=1^+(1^{+-})$ charm and bottom sectors in the framework of a constituent quark model which satisfactorily describes a wide range of properties of (non-)conventional hadrons containing heavy quarks. All the relevant channels are included for each sector, i.e.: The $D^{(\ast)}\bar D^{\ast}+h.c.$, $\pi J/\psi$ and $\rho\eta_c$ channels for the $Z_c$ and $B^{(\ast)}B^{\ast}$ and $\Upsilon(nS)\pi$ ($n=1,2,3$) channels for the $Z_b$ analysis. Possible structures of these resonances will be discussed.


Introduction
The search of exotic structures in the heavy meson spectra, beyond the simplest qq structure, is relatively recent. In 2003, the first and most iconic resonance, the so-called X (3872), was spotted by Belle [1] and promptly confirmed by other B-factories and accelerator facilities such as BaBar, CDF and D0 Collaborations. At the same time, BaBar and CLEO Collaborations reported the discovery of the D * s0 (2317) and D s1 (2460) [2,3], two puzzling structures in the heavy-light mesons spectrum. Even though the properties of such resonances where compatible with a qq quark content, they were difficult to accommodate in the naive quark model due to their unexpected masses and decay properties, which pointed to a non-negligible contribution of higher Fock-state components.
The clearest evidence of exotic structures appeared in 2011, when the Belle Collaboration [4,5] announced the observation of two meson-like structures in the bottom sector with forbidden quantum numbers for a bb pair. The so-called Z b (10610) and Z b (10650) where charged structures close to the BB * and B * B * thresholds, respectively, spotted in the Υ (5S) → π + π − Υ (nS) reaction. Few years later, their charmonium partners arrived. The BESIII and Belle Collaborations discovered the Z c (3900) [6,7], close to the DD * threshold, in the π + π − J/ψ invariant mass spectrum of the e + e − → π + π − J/ψ process at s = 4.26 GeV. The spin and parity of the charged Z c (3900) was set to J P = 1 + [8], with a mass M = (3881.2 ± 4.2 ± 52.7) and width Γ = (51.8 ± 4.6 ± 36.0) MeV. Additionally, the Z c (4020) resonance, close to the D * D * threshold, was discovered soon after at BESIII in the e + e − → π + π − h c channel with a mass of M = (4022.9 ± 0.8 ± 2.7) MeV/c 2 and a width of Γ = (7.9 ± 2.7 ± 2.6) MeV [9].
All the previous Z b and Z c states are charged, so they cannot be described as pure qq states. Their closeness to B ( * )B * and D ( * )D * thresholds indicate a dominant meson-meson component in their wave functions. These features have been widely explored in different theoretical scenarios, such as hadron molecules [10][11][12][13], tetraquark structures [14][15][16][17] or simple kinematic effects linked to the opening of meson-meson thresholds [18,19].
In this work we explore the I G (J P C ) = 1 + (1 +− ) charm and bottom sector in a coupled-channels scheme, including the closest meson-meson thresholds. The meson-meson interaction is described in the framework of a constituent quark model 1 successfully employed to explain the meson and baryon phenomenology from the light to the heavy quark sector (see for example Ref. [20,[22][23][24]. Moreover, the D ( * )D * residual interaction deduced from the model has been satisfactorily used to describe meson-meson [25][26][27] molecular states.

Theoretical Formalism
Invariance under chiral rotations is a symmetry of the Quantum Chromodynamics (QCD) Lagrangian with massless light quarks. However, this symmetry is not satisfied in nature, pointing to a spontaneous breaking within QCD. This effect has several interesting consequences, one of those being the emergence of a momentum-dependent constituent quark mass, M = M (q 2 ) and M (q 2 → 0) = m q , with m q the current quark mass, and the appearance of Goldstone bosons which mediate the interaction among light quarks.
The latter phenomenology can be synthesized in a constituent quark model (CQM), which is described by the following low-energy Lagrangian [28] L where U γ 5 = e iλ a φ a γ 5 / f π is the Goldstone-boson fields matrix. This matrix can be expanded as so the first term can be identified as the constituent quark mass, the second one describes the pseudoscalar meson exchange interaction among quarks and the third term, whose main contribution is the two-pion exchange, can be coded by means of a scalar-meson exchange potential. The model is completed with QCD non-perturbative effects such as the confinement interaction, implemented phenomenologically so colored hadrons are prohibited. Within our CQM, the confinement is modelled with a linear-rising potential, due to multi-gluon exchanges among quarks, which is screened at large inter-quark distances due to sea quarks [29]: Here, a c and µ c are model parameters. From the latter equation it is easy to see that the potential is linear at short inter-quark distances with an effective confinement strength Beyond the non-perturbative energy scale, the dynamics of the quarkonium is expected to be dominated by QCD perturbative effects. That is taken into account by means of the one-gluon exchange potential, derived from the following Lagrangian Here, α s is the strong coupling constant, λ a are the SU(3) colour matrices and G µ a is the gluon field. To consistently treat the light-and heavy-quark sectors we employ a gluon coupling constant that scales with the reduced mass of the interacting quarks. Explicit expressions, model parameters and a detailed physical background of the constituent quark model used in this work can be found in, e.g., Refs. [20,23] These CQM describes the interaction among constituent quarks. When dealing with meson (hadron) degrees of freedom, the Resonating Group Method [30] can be employed to obtain the interaction at meson level from the microscopic interaction at quark level. In RGM, mesons are considered as quark-antiquark clusters and effective cluster-cluster interaction emerges from the underlying qq dynamics.
The meson eigenstates φ C ( p C ) of a general meson C, with p C the relative momentum between the quark and antiquark of the meson C, are calculated by means of the two-body Schrödinger equation using the Gaussian Expansion Method [31]. We use Gaussian trial functions with ranges given by a geometrical progression [31], which optimizes the method and reduces the number of free parameters. This choice produces a dense distribution at short distances, enabling a better description of the dynamics mediated by short range potentials.
The orbital wave function of a system composed of two mesons A and B with distinguishable quarks can be, then, written as 2 where χ α ( P) is the relative wave function between the two clusters, and α labels the set of quantum numbers needed to uniquely define a certain partial wave.
Such relative wave function can be obtained as the solution of the projected Schrödinger equation: where E is the energy of the system. As we can see, two type of diagrams appear: one which does not consider quark exchanges between clusters, called direct potential; and one involving quark exchanges, dubbed rearrangement potential. The direct potential RGM V αα D ( P , P) of a reaction AB → C D can be written as where {i, j} runs over the constituents of the involved mesons. The quark rearrangement potential RGM V αα R ( P , P) represents a natural way to connect meson-meson channels with different quark content, such as πJ/ψ and DD * . It can be calculated with where P mn is the operator that exchanges quarks between clusters. The solution of the RGM coupled-channels equation is performed deriving from Eq. (6) a set of coupled Lippmann-Schwinger equations of the form Here, 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 non-relativistic case as: where µ α is the AB-system reduced mass corresponding to the channel α, and ∆M α is the difference between the threshold of the AB system and the reference one. The coupled-channels Lippmann-Schwinger equation is solved via the matrix-inversion method proposed in Ref. [32], generalized in order to include channels with different thresholds. From the full T -matrix, it is direct to extract the on-shell part, directly related to the scattering matrix as in non-relativistic kinematics, with k α the on-shell momentum for channel α.
The final goal of the study is the analysis of the existence of states above and below thresholds within the same formalism. For that purpose, all the potentials and kernels are analytically continued for complex momenta so the poles of the T -matrix in any possible Riemann sheet can be detected.

Z c (3900) and Z c (4020) structures
The aim of the present study is to use the constituent quark model of Ref. [20] as a basis to perform a coupled-channels calculation of the I G (J P C ) = 1 + (1 +− ) sector with hidden charm and bottom. For the charmonium sector, we include the closest thresholds to the Z c (3900) and Z c (4020) experimental masses: πJ/ψ (3234.19 MeV/c 2 ), ρη c (3755.79 MeV/c 2 ), DD * (3875.85 MeV/c 2 ), D * D * (4017.24 MeV/c 2 ), where the threshold masses are shown in parenthesis. The h c π channel is not considered as it couples weakly to the D ( * )D * channels. In fact, the only contribution of this channel would come from the 3 D 1 component of the internal wave function of the D * meson, which in our model is ∼ 0.03% and it is neglected in the present calculation. For the 3 S 1 component of the D * , the coupling to h c π is exactly zero. Similarly, we do not include other nearby channels such as the χ cJ ρ ones, whose couplings are found to be almost three orders of magnitude smaller than those of the J/ψπ and η c ρ channels.
The invariant mass distribution of the DD * , πJ/ψ and D * D * channels predicted by our model are shown in Figs. 1 and 2, following the procedure of Ref. [33]. Normalization factors N AB and production amplitudes A AB , that appear in the line shapes, are obtained from a fit performed with the experimental results of Refs. [6,8,34], that is, experimental data for DD * and πJ/ψ channels. The amplitudes A AB obtained from such fit are, then, employed for the D * D * channel in order to obtain the global normalization N D * D * from experimental data of Ref [35]. In order to describe the experimental measurement, the theoretical line shapes have been convoluted with the detector resolution.
In order to optimize the fit, only DD * experimental data up to 3.92 GeV was considered. For larger energies the background is expected to dominate [36]. The πJ/ψ data of Ref. [8] was fitted from 3.85 GeV on. In the theoretical line shape of the DD * channel one can clearly see two enhancements related with the opening of the DD * and D * D * thresholds, which can be associated with the Z c (3900) and the Z c (4020). In the πJ/ψ case, only one enhancement appear around 3.87 GeV while the opening of the D * D * channel appears as a slight step down in the number of events.
Apart from the unknown π + A + B vertex details, encoded in the fitted normalization and amplitude factors, our CQM is able to reproduce the experimental data with no fine-tuning of the parameters. Attending to the nature of the Z c (3900) and the Z c (4020) states, we have examined the pole structure of the S-matrix for different coupled-channels calculations. Our results are shown in Table 1. For the Z c (3900), even for a one-channel DD * calculation, the S-matrix shows a virtual pole below threshold, which is maintained when further channels are added. In the complete calculation, the Z c (3900) is associated with a pole located in the imaginary axis of the second Riemann sheet below the DD * threshold and, thus, it is a virtual state. The situation is similar in the case of the Z c (4020), which is interpreted as a virtual state located below the D * D * threshold. Further details of the calculation can be found in Ref. [33].

Z b (10610) and Z b (10650) structures
For the bottomonium sector we follow the same approach as above, i.e., we perform a coupledchannels calculation of the I G (J P C = 1 + (1 +− ) sector including the following channels: BB * produce a real bound state below the BB * threshold and a resonance close to the B * B * threshold, as can be seen in Table 2. Due to Heavy Flavour Symmetry, the B ( * )B * and D ( * )D * interactions are practically identical. However, the kinetic energy of the involved channels is reduced in the bottom sector due to the larger mass of the b-quark mass, favouring the creation of bound states. In order to see if such predicted poles describe the experimental situation, analogously as for the Z c 's states, we calculated the predicted line shapes in each channel between 10.55 and 10.75 GeV. We consider they are coming from a π + A + B vertex with a center of mass energy of s = 10.865 GeV, according to Refs. [4,5,37]. The procedure to obtain the line shapes is the same as for the Z c 's (see Ref. [33]). The results can be seen in Figs. 3 for the BB * , B * B * , Υ (1S)π and Υ (2S)π. The description of the Z b (10610) peak is in nice agreement with the experimental data for the BB * and Υ (1S)π channels. The Z b (10650) peak is also properly described for the B * B * channel, but somehow it is missing in the BB * and Υ (1S)π line shapes. Such discrepancy could be caused by the simple ansatz considered for the π + A + B vertex or perhaps it points to a small  Table 1: The S-matrix pole positions, in MeV/c 2 , for different coupled-channels calculations. The included channels for each case are shown in the first column. Poles are given in the second and fourth columns by the value of the complex energy in a specific Riemann sheet (RS). The RS columns indicate whether the pole has been found in the first (F) or second (S) Riemann sheet of a given channel. Each channel in the coupled-channels calculation is represented as an array's element, ordered with increasing energy.  non-diagonal coupling of the B * B * channel with the rest of the thresholds within our model which should be improved. Both peaks emerge in the Υ (2S)π line shape, but their strength is not high enough to describe the experimental situation. Further research is ongoing in order to clarify the inner structure of the Z b (10610) and Z b (10650).

Conclusion
The Z b 's and Z c 's are very peculiar structures, different from other molecular states of the bottomonium and charmonium spectrum. In order to clarify their inner structure, we have performed, within the framework of a constituent quark model, a coupled-channels calculation of the I G (J P C ) = 1 + (1 +− ) sector around the energies of the Z c (3900) and Z c (4020), for the charmonium sector, and the Z b (10610) and Z b (10650), for the bottomonium case, including the most relevant thresholds. The line shapes of the D ( * )D * , πJ/ψ, B ( * )B * invariant mass distributions are well reproduced without any fine-tuning of the model parameters. Some channels such as the Υ (nS)π are not well understood yet and will require further investigation. The analysis of the S-matrix poles allows us to conclude the following. For the Z c 's the point-wise behavior of the line shapes is due to the presence of two virtual states that can be seen as D ( * )D * threshold bumps, and have an overall good agreement with the location of the Z c (3900) and Z c (4020) signals. These results confirm the conclusion of the lattice QCD calculation of Ref. [38]. For the Z b 's, a real bound state, below the BB * threshold, and a resonance, below the B * B * threshold, emerge from the coupled-channels interactions, matching almost all the properties of the Z b (10610) and Z b (10650) structures, respectively.