The $Z_c$ structures in a coupled-channels model

The $Z_{c}(3900)^\pm/Z_c(3885)^\pm$ and $Z_{c}(4020)^\pm$ are two charmonium-like structures discovered in the $\pi J/\psi$ and $D^\ast\bar D^{(\ast)}+h.c.$ invariant mass spectra. Their nature is puzzling due to their charge, which forces its minimal quark content to be $c\bar c u\bar d$ ($c\bar c 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 J/\psi$ and the closeness of their mass to $D^\ast\bar D^{(\ast)}$-thresholds stimulates both a molecular interpretation or a coupled-channels threshold effect. In this work we perform a coupled-channels calculation of the $I^G(J^{PC})=1^+(1^{+-})$ sector including $D^{(\ast)}\bar D^{\ast}+h.c.$, $\pi J/\psi$ and $\rho\eta_c$ channels in the framework of a constituent quark model which satisfactorily describes a wide range of properties of (non-)conventional hadrons containing heavy quarks. The meson-meson interactions are dominated by the non-diagonal $\pi J/\psi-D^\ast\bar D^{(\ast)}$ and $\rho\eta_c-D^\ast\bar D^{(\ast)}$ couplings which indicates that the $Z_{c}(3900)^\pm/Z_c(3885)^\pm$ and $Z_{c}(4020)^\pm$ are unusual structures. The study of the analytic structure of the $S$-matrix allows us to conclude that the point-wise behavior of the line shapes in the $\pi J/\psi$ and $D\bar D^*$ invariant mass distributions is due to the presence of two virtual states that produce the $Z_c$ peaks.


I. INTRODUCTION
Since early 2000s, signals of non-conventional meson structures have appeared in the so-called B-factories and other accelerator facilities. In 2003, the X(3872) was discovered by Belle [1], promptly confirmed by BaBar [2], CDF [3] and D0 [4] Collaborations, whereas BaBar and CLEO Collaborations reported the observation of two puzzling heavy-light mesons: D * s0 (2317) and D s1 (2460) [5,6]. The common characteristic of all these states is that, although their quantum numbers are compatible with naive qq structures, their masses and decay properties point to more complex structures involving higher Fock-state components.
However, it was not until 2011 when undeniable evidences of exotic mesons with forbidden quantum numbers for a quark-antiquark pair were observed. These indications were found by the Belle Collaboration [7,8] in the bottom sector with the discovery of the Z b (10610) and Z b (10650) charged structures in the Υ(5S) → π + π − Υ(nS) reaction. As the reader can guess, they are close to the BB * and B * B * thresholds, respectively.
One of the partners of such structures in the charmonium sector arrived few years later, in 2013, when the BESIII and Belle Collaborations claimed the observation of a charged structure, dubbed Z c (3900) [9,10], in the π + π − J/ψ invariant mass spectrum of the e + e − → π + π − J/ψ process at √ s = 4.26 GeV. It was further confirmed in the same channel but at √ s = 4.17 GeV in Ref. [11] by a reanalysis of CLEO-c data. Its neutral partner was also reported in Refs. [11,12].
The Z c (3900) mass from the πJ/ψ invariant mass spectrum measurement, M = (3899.0 ± 3.6 ± 4.9) MeV/c 2 [9], is 24 MeV/c 2 above the DD * mass threshold and thus it is natural to observe the Z c (3900) in the (DD * ) I=1 channel. This was confirmed by the BESIII Collaboration [13] which reported a strong threshold enhancement in the DD * invariant mass distribution of the e + e − → π ± (DD * ) ∓ process at √ s = 4.26 GeV, referred as Z c (3885). Fitted with a mass-dependentwidth Breit-Wigner line shape, its mass and width were determined to be M = (3883.9 ± 1.5 ± 4.2) MeV/c 2 and Γ = (24.8 ± 3.3 ± 11.0) MeV, the data being compatible with the one expected for J P = 1 + quantum numbers. A few years later, such determination was updated to M = (3881.7 ± 1.6 ± 6.9) MeV/c 2 and Γ = (26.6 ± 2.0 ± 2.1) MeV, with measurements at √ s = 4.26 GeV and √ s = 4.23 GeV [14]. Recently, the spin and parity of the Z c (3900) was analyzed in Ref. [15] with the result that the J P = 1 + assignement is favoured by the data; moreover, its mass was reduced to M = (3881.2 ± 4.2 ± 52.7) and the new width update to Γ = (51.8 ± 4.6 ± 36.0) MeV, using a Flatté-like expresion. This new determination supports the idea that the Z c (3900)/Z c (3885) signals correspond to an unique state.
Soon after all these experimental activities, the BESIII Collaboration reported the discovery of another charged state, the Z c (4020) resonance, 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 [16]. Later on, a similar structure with a mass of M = (4026.3 ± 2.6 ± 3.7) MeV/c 2 and a width of Γ = (24.8 ± 5.6 ± 5.7) MeV was observed by the BESIII Collaboration in the process e + e − → π ± (D * D * ) ∓ at √ s = 4.26 GeV [17]. It is worth also to note that the neutral partner of the Z c (4020) was reported by the BESIII Collaboration in Ref. [18].
The fact that most of these structures have been found at √ s = 4.26 GeV led several authors to think that the Y (4260) resonance may be responsible for the production of the Z c (3900) and Z c (4020) structures due to their unconventional properties. However, the detection of the Z c (3900) at the peak of the ψ(4160) state by the CLEOc Collaboration discarded any effect of the Y (4260) resonance [11].
The Z c (3900) and Z c (4020) states are charged and close to DD * and D * D * thresholds, respectively. Therefore, they cannot be described as pure qq states and it is expected that D ( * )D * components are significant in their wave functions. This fact can be implemented in different theoretical scenarios, ranging from hadron molecules [19][20][21][22], to tetraquark structures [23][24][25][26][27][28] or simple kinematic effects linked to the opening of mesonmeson thresholds [29,30].
At this point, it is worth noticing that not all enhancements of the cross section correspond to a resonance. As discussed by R. G. Newton [31], the vicinity of the threshold of a new opening channel causes the cross section to show an anomalous behavior in form of a cusp. From a mathematical point of view, this behavior is produced by the branch-point of the S-matrix at threshold. Resonances are also enhancements of the cross section, but corresponding to poles of the S-matrix on the unphysical Riemann sheet close to the physical one. In this case, the cross section can be parametrized with a Breit-Wigner form with two parameters: the energy (mass) and the width. Furthermore, besides bound states which correspond to poles of the S-matrix on the real energy axis of the physical Riemann sheet below the lowest threshold, the S-matrix may present poles on the real energy axis but in the unphysical Riemann sheet. If these virtual states have sufficiently low energy, then their presence modifies appreciably the energy dependence of the scattering cross section which can be anomalously large. Therefore, to correctly interpret the point-wise behavior of the line shapes of a particular process one must look carefully on the analytical properties of the S-matrix.
In a couple of papers, E. S. Swanson [29,30] claimed that the Z c and Z b resonances could arise from kinematic threshold effects, meaning that the DD * interaction, in case of Z c (3900), is strong enough to cause an enhancement just above the S-wave threshold. The fact that the masses of such resonances are slightly above their respective thresholds seems to favor this possibility. Swanson used a simplified model to describe Belle and BESIII data as kinematic cusps, successfully capturing the features of all the data and indicating that there is no evidence for strong DD * or D * D * rescattering in the (I = 1) J P C = 1 +− channel. Hence, he concluded that isovector rescattering is not sufficiently attractive to generate dynamical bound states and so exotic resonances are not required to explain the data.
Kinematic cusp interpretations were subsequently challenged by F.-K. Guo et al. [32]. Using a nonrelativistic effective field theory, the authors argued that a cusp always exists due to the opening of an S-wave threshold but, in order to reproduce the narrow and pronounced peak observed by experimentalists, a nonperturbative interaction amongst the heavy mesons is necessary and thus a nearby pole in the S-matrix must appear. This motivated a detailed analysis of the experimental data using reaction theory and considering different scenarios for the inner structure of the Z c (3900) [33]. The authors claimed that a cusp is not strong enough to reproduce the peak while a molecular or virtual state can describe equally well the experimental data.
Further calculations with different effective field theories have been performed. F. Aceti et al. [34], using the local hidden gauge approach, studied the DD * interaction in a coupled channel Bethe-Salpeter calculation. They found in the I = 1 channel a barely DD * bound state decaying into ρη c and πJ/ψ channels with a mass of 3869 − 3875 MeV. Additionally, M. Albaladejo et al. [35] performed a simultaneous study of the invariant mass distribution for the πJ/ψ and DD * channels and found different interpretations when exploring with an energy dependent or independent DD * S-wave interaction. In the first case, the Z c enhancement was originated from a resonance with a mass around the DD * threshold. In the second one, the Z c signal was produced by a virtual state which must have a hadronic molecular nature. In any case, a DD * bound state was not favored. It is worth remarking that, following these authors, the tetraquark scenario would be discarded if the S-matrix pole is at the unphysical sheet (virtual state), as it would originate from DD * interaction solely. A similar conclusion can be found in Ref. [36], [21] where the authors look for the Z c (3900) and Z c (4020) state in the πJ/ψ and DD * and the D * D * invariant mass spectra respectively. The calculation was done in a quasi-potential Bethe-Salpeter equation approach at the hadronic level. They found that the πJ/ψ-DD * interaction produces a virtual state at a mass around 3870 MeV whereas they found abound state for the Z c (4020) with reasonable values of the parameters.
Many papers have been released in the last years reporting exploratory LQCD studies of the Z c structures [37][38][39][40][41][42]. S. Prelovsek et al. [37,38] looked for the Z c (3900) + state in the energy region below 4 GeV, finding only two-particle πJ/ψ and DD * scattering states and no signal of the Z c (3900) + . In a similar search Y. Chen et al. [39] found that the DD * interaction is weakly repulsive and, therefore, the results do not support the possibility of a shallow bound state at least at the pion mass values studied .A different approach is taken by Y. Ikeda et al. [42]. They performed a coupledchannels calculation taking into account the πJ/ψ, ρη c and DD * channels and finding that the interactions between them are dominated by off-diagonal couplings which may indicate that the Z c (3900) + can be explained as a threshold cusp. There are technicalities involving LQCD computations, such as large pion masses, small volumes and a set of interpolators not large enough for having overlap with the physical state, which still prevent to make a definitive statement. Nonetheless, seems that the available LQCD calculations are robust enough to discard the bound state option for the Z + c states. From a molecular point of view, the fact that the Z c structures have I = 1 makes that the D ( * )D( * ) interaction weak because in these channels the onepion-interacion is mainly repulsive. This leaves the coupled channel calculations as the most promise option to produce some resonance, bound or virtual state if any In this work we analyze the I G (J P C ) = 1 + (1 +− ) 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 [45][46][47][48][49][50][51][52][53][54]. Moreover, the D ( * )D * interaction deduced from the model has been satisfactorily used to describe meson-meson [55][56][57] and meson-baryon [58][59][60][61] molecular states such as the X(3872) as a DD * molecule coupled to cc(n 3 P 1 ) states [62,63].
The structure of the present manuscript is organized in the following way. In Sec. II the theoretical framework is briefly presented and discussed. Section III is devoted to the analysis and discussion on the obtained results. We summarize and give some conclusions in Sec. IV.

II. THEORETICAL FORMALISM
A. Constituent quark model The Lagrangian of Quantum Chromodynamics (QCD) with massless light quarks is invariant under chiral rotations. This symmetry does not appear in Nature indicating that it is spontaneously broken in QCD. Among other consequences, a constituent mass which depends on the quark momentum, M = M (q 2 ) and M (q 2 = ∞) = m q , is developed and Goldstone-boson exchange interactions appear between the light quarks.
Our constituent quark model (CQM) tries to mimic the previous phenomena based on the following effective 1 The interested reader is referred to Refs. [43,44] for detailed reviews about the naive quark model in which this work is based.
Lagrangian at low-energy [64] L being U γ5 = e iλaφ a γ5/fπ the matrix of Goldstone-boson fields that can be expanded as The constituent quark mass is obtained from the first term, the second one describes the pseudoscalar meson exchange interaction among quarks and the main contribution of the third term comes from the two-pion exchange which is modeled by means of a scalar-meson exchange potential. Another nonperturbative effect is the confining interaction which is implemented phenomenologically in order to avoid colored hadrons. In our CQM, the confinement is represented as a linear potential, due to multi-gluon exchanges between quarks, that is screened at large interquark distances, as a consequence of sea quarks [65]: Here, a c and µ c are model parameters. One can see that the potential is linear at short inter-quark distances with an effective confinement strength σ = −a c µ c ( λ c i · λ c j ), while it becomes constant at large distances.
Beyond the nonperturbative energy scale one expects that the dynamics of the bound-state system is governed by QCD perturbative effects. We take it into account through the one-gluon exchange potential derived from the following vertex Lagrangian where α s is the strong coupling constant, λ a are the SU (3) colour matrices and G µ a is the gluon field. The strong coupling constant, α s , has a scale dependence which allows a consistent description of light, strange and heavy mesons. Its explicit expression can be found in, e.g., Ref. [45].
A detailed physical background of the quark model can be found in Refs. [45,66]. 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 flavor-blind.

B. Resonating Group Method and
Lippmann-Schwinger equation The aforementioned CQM specifies the microscopic interaction among constituent quarks. In order to describe the interaction at the meson level we employ the Resonating Group Method [67], where mesons are considered as quark-antiquark clusters and an effective cluster-cluster interaction emerges from the underlying qq dynamics.
We assume that the wave function of a system composed of two mesons A and B with distinguishable quarks can be written as 2 where φ C ( p C ) is the wave function of a general meson C with p C the relative momentum between the quark and antiquark of the meson C. The wave function which takes into account the relative motion of the two mesons is χ α ( P ). The projected Schrödinger equation for the relative wave function can be written as follows: where E is the energy of the system. The direct potential RGM V αα D ( P , P i ) can be written as The quark rearrangement potential RGM V αα R ( P , P i ) represents a natural way to connect meson-meson channels with different quark content, such as πJ/ψ and DD * , and it is given by where P mn is the operator that exchanges quarks between clusters. The meson eigenstates φ C ( p C ) are calculated by means of the two-body Schrödinger equation, using the Gaussian Expansion Method [68]. This method 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, we use Gaussian trial functions whose ranges are given by a geometrical progression [68]. This choice produces a dense distribution at short distances enabling better description of the dynamics mediated by short range potentials.
The solution of the coupled-channels RGM equations is performed deriving from Eq. (6) a set of coupled Lippmann-Schwinger equations of the 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: Here, µ α is the reduced mass of the AB system corresponding to the channel α, and ∆M α is the difference between the threshold of the AB system and the one we take as a reference.
We solve the coupled-channels Lippmann-Schwinger equation using the matrix-inversion method proposed in Ref. [69], generalized 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 nonrelativistic kinematics): with k α the on-shell momentum for channel α.
Our aim is to explore the existence of states above and below thresholds within the same formalism. Thus, we have to analytically continue all the potentials and kernels for complex momenta in order to find the poles of the T -matrix in any possible Riemann sheet.

C. Production line shapes of the Zc
Most of the experimental data of the Z c (3900) and Z c (4020) were taken at √ s = 4.26 GeV, which corresponds to the energy of the resonance Y (4260), though there are several measurements at different energies such as √ s = 4.17 GeV [11]. To simplify the theoretical description of the Z c generating process and to avoid possible complexities of considering the internal structure of the Y (4260) resonance, we adopt a phenomenological vertex which creates the π + A + B triplet, where A + B is the final two-body state of the coupled-channels calculation, i.e. A + B = {πJ/ψ, ρη c , DD * , D * D * }. The three-body decay π + A + B of a resonance in its rest frame can be written as where M is the production amplitude and dΦ is the three-body phase-space given by Therefore, we can express dΓ in terms of the invariant mass spectrum of the two-meson channel, m AB , as where (k AB , Ω AB ) is the on-shell momentum of the AB pair and Ω πZc is the solid angle of the π in the center-ofmass rest frame of π + A + B at energy √ s. The on-shell momenta are given by where λ(M, m 1 , , with m ± = m 1 ± m 2 . Integrating the angles, we have with β the quantum numbers of the channel AB. The Lorentz-invariant production amplitude, M, describes the Z c → AB production. It is diagramatically shown in Fig. 1 and can be written as (18) where β ( ) denotes a given AB channel in the coupledchannels calculation.
We only consider production through the AB in S-wave (as the vertex does not have momentum dependence) and we add the amplitudes A β to take into account different production probabilities for each channel AB.
Experimentalists measure events in the σ(e + e − → π ± Z ∓ c ) × B(Z ∓ c → AB) process and thus, to describe the data, we need to add a normalization factor to translate the decay rate into events: This normalization encodes other relevant process details such as the value of σ(e + e − → πZ c ) if the Z c In order to fit {A AB , N AB } for each channel, we minimize a global χ 2 function using all the available experimental data for channels πJ/ψ, DD * and D * D * : and the uncertainty of the parameters {A, N } are estimated from the second derivative of χ 2 at its minimum. That is, for example, around the best-fit value N BF , the χ 2 function can be approximated by with σ N the error of the best-fit value N BF .

III. RESULTS
As already mentioned above, we perform a coupledchannels calculation within the framework of the constituent quark model of Ref. [45] for the I G (J P C ) = 1 + (1 +− ) sector, including the closest thresholds to the experimental masses of the Z c (3900) and Z c (4020) hadrons, that is: π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.
Our results for the invariant mass distribution of the DD * , πJ/ψ and D * D * channels are shown in, respectively, the upper and lower panels of Fig. 2 and in Fig. 3. Note that all channels mentioned in the above paragraph are included in the calculation of the theoretical line shapes. To translate the decay rates into events, we use a normalization factor N AB fitted for each channel according to the expression of Eq. (20). The amplitudes A AB at Eq. (18) are the same for all the line shapes. The shaded area around the theoretical curve shows the statistical 68%-confident level (CL) of the fit, obtained by propagating the errors of the fitted parameters by means of the covariance matrix. The χ 2 fit is performed with the experimental results of Refs. [9,14,15], 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 [17]. Line shapes for DD * (upper pannel) and πJ/ψ (lower pannel) at √ s = 4.26 GeV. Experimental data are from Ref. [14,15], respectively. The theoretical line shapes have been convoluted with the experimental resolution. The line-shape's 68% uncertainty is shown as a shadowed area.
shown in Table I, the result on the χ 2 /d.o.f. is also collected therein. In order to describe the experimental measurement, the theoretical line shapes have been convoluted with the detector resolution.
For the fit, the DD * data was taken from the m DD * threshold up to 3.92 GeV. This is the region where the signal of the Z c (3900) dominates over the background [13], from there on the data set is too noisy. The πJ/ψ data of Ref. [15] 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.
Besides the unavoidable normalization factors and amplitudes related to the phenomenological π + A + B vertex, our formalism is able to reproduce the experimental data without fine-tuning any model parameter. The meson-meson interaction driven by the quark rearrangement process, which connect the D ( * )D * channels with the πJ/ψ and ρη c ones, is of the same order of magnitude than the direct interaction, D ( * )D * − D ( * )D * . This points to an important role of the πJ/ψ(ρη c ) − D * D( * ) mixings in this sector, which could be an indication of the kinematic nature of the Z c states as some authors claimed recently [36,42].
To disentangle the contribution of the different channels to the line shape, we compare in Fig. 4 the line shape for the full calculation with those taking into account partial combinations of the considered channels. When the DD * scattering channel is considered alone, a small enhancement basically due to the pion tensor interaction between the S and D waves appears but the generated peak is too wide. The inclusion of the ρη c channel narrows that peak, making it more compatible with the experimental situation. On the other hand, adding to the DD * channel the D * D * one generates a second structure at its threshold opening associated to the Z c (4020) peak. One can see that the second enhancement is much higher  Table I. See text for discussion.
than the experimental data if we only include D ( * )D * channels. The line shape moves closer to the experimental situation when the ρη c and πJ/ψ channels are considered in the calculation, showing that they play an important role in building the observed enhancements.
To deepen into the nature of the Z c (3900) and the Z c (4020), we have examined the analytic structure on the complex energy plane of the S-matrix for the different coupled-channels calculations. Our results are shown in Table II. For the Z c (3900), one can see that even for a one-channel DD * calculation the S-matrix shows a pole below threshold in the second Riemann sheet. When the ρη c is included a pole corresponding with a virtual state appears. The inclusion of the rest of the channels does not change drastically the pole position. However, it is necessary to reproduce the experimental line shape. In the complete calculation, the Z c (3900) is associated with a pole located in the imaginary axis and 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.
It is worth noticing that the real part of the Z c (3900) mass pole, shown in Table II, is always below the DD * threshold and, hence, it does not seem to be compatible with the mass and width estimations from experimental measurements [9,14,15]. Obviously, in our calculation, the Z c (3900) and Z c (4020) are not resonances but virtual state and, therefore, a directly comparison of the complex energy of the pole with the Breit-Wigner parametrization of a resonance would be misleading. The only way to connect our results with the experimental data is through the description of the line shapes.

IV. SUMMARY
The Z c (3900) and Z c (4020) are very peculiar structures different from other molecular resonances of the charmonium spectrum. As they have I = 1, 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 channels responsible for the structures appeared in the line shapes. Therefore, they should only appear in calculations involving coupled channels.
To show up this mechanism, 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), including the most relevant thresholds. The lines shapes of the DD * , πJ/ψ and D * D * invariant mass distributions are well reproduced without any fine-tuning of the model parameters. The analysis of the S-matrix poles allows us to conclude that the point-wise behavior of the line shapes is due to the presence of two virtual states which produces the Z c (3900) and Z c (4020) signals. These results confirm the conclusion of the lattice QCD calculation of Ref. [42].  II. 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.