Eur. Phys. J. C (2019) 79:78 https://doi.org/10.1140/epjc/s10052-019-6552-7 Regular Article - Theoretical Physics The Zc structures in a coupled-channels model Pablo G. Ortega1,a , Jorge Segovia2,b , David R. Entem1,c , Francisco Fernández1,d 1 Grupo de Física Nuclear and Instituto Universitario de Física Fundamental y Matemáticas (IUFFyM), Universidad de Salamanca, 37008 Salamanca, Spain 2 Grup de Física Teòrica, Dept. Física and IFAE-BIST, Universitat Autònoma de Barcelona, 08193 Bellaterra, Barcelona, Spain Received: 3 August 2018 / Accepted: 31 December 2018 / Published online: 28 January 2019 © The Author(s) 2019 Abstract The Z c (3900)± /Z c (3885)± and Z c (4020)± are characteristic of all these states is that, although their quan- two charmonium-like structures discovered in the π J/ψ and tum numbers are compatible with naive q q̄ structures, their D ∗ D̄ (∗) +h.c. invariant mass spectra. Their nature is puzzling masses and decay properties point to more complex struc- due to their charge, which forces its minimal quark content to tures involving higher Fock-state components. be cc̄u d̄ (cc̄d ū). Thus, it is necessary to explore four-quark However, it was not until 2011 when undeniable evidences systems in order to understand their inner structure. Addi- of exotic mesons with forbidden quantum numbers for a tionally, their strong coupling to channels such as π J/ψ quark–antiquark pair were observed. These indications were and the closeness of their mass to D ∗ D̄ (∗) -thresholds stim- found by the Belle Collaboration [7,8] in the bottom sector ulates both a molecular interpretation or a coupled-channels with the discovery of the Z b (10610) and Z b (10650) charged threshold effect. In this work we perform a coupled-channels structures in the ϒ(5S) → π + π − ϒ(nS) reaction. As the calculation of the I G (J PC ) = 1+ (1+− ) sector including reader can guess, they are close to the B B̄ ∗ and B ∗ B̄ ∗ thresh- D (∗) D̄ ∗ +h.c., π J/ψ and ρηc channels in the framework of a olds, respectively. constituent quark model which satisfactorily describes a wide One of the partners of such structures in the charmonium range of properties of (non-)conventional hadrons containing sector arrived few years later, in 2013, when the BESIII and heavy quarks. The meson–meson interactions are dominated Belle Collaborations claimed the observation of a charged by the non-diagonal π J/ψ − D ∗ D̄ (∗) and ρηc − D ∗ D̄ (∗) structure, dubbed Z c (3900) [9,10], in the π + π − J/ψ invari- couplings which indicates that the Z c (3900)± /Z c (3885)± ant mass spectrum of the e+ e− → π + π − J/ψ process at √ and Z c (4020)± are unusual structures. The study of the ana- s = 4.26 GeV. It was further confirmed in the same chan- √ lytic structure of the S-matrix allows us to conclude that the nel but at s = 4.17 GeV in Ref. [11] by a reanalysis of point-wise behavior of the line shapes in the π J/ψ and D D̄ ∗ CLEO-c data. Its neutral partner was also reported in Refs. invariant mass distributions is due to the presence of two vir- [11,12]. tual states that produce the Z c peaks. The Z c (3900) mass from the π J/ψ invariant mass spec- trum measurement, M = (3899.0 ± 3.6 ± 4.9) MeV/c2 [9], is 24 MeV/c2 above the D D̄ ∗ mass threshold and thus it 1 Introduction is natural to observe the Z c (3900) in the (D D̄ ∗ ) I =1 chan- nel. This was confirmed by the BESIII Collaboration [13] Since early 2000s, signals of non-conventional meson struc- which reported a strong threshold enhancement in the D D̄ ∗ tures have appeared in the so-called B-factories and other invariant mass distribution of the e+ e− → π ± (D D̄ ∗ )∓ √ accelerator facilities. In 2003, the X (3872) was discovered process at s = 4.26 GeV, referred as Z c (3885). Fitted by Belle [1], promptly confirmed by BaBar [2], CDF [3] with a mass-dependent-width Breit-Wigner line shape, its and D0 [4] Collaborations, whereas BaBar and CLEO Col- mass and width were determined to be M = (3883.9 ± laborations reported the observation of two puzzling heavy- 1.5 ± 4.2) MeV/c2 and  = (24.8 ± 3.3 ± 11.0) light mesons: Ds0∗ (2317) and D (2460) [5,6]. The common MeV, the data being compatible with the one expected s1 for J P = 1+ quantum numbers. A few years later, a e-mail: pgortega@usal.es such determination was updated to M = (3881.7 ± 1.6 b e-mail: jsegovia@ifae.es ± 6.9) MeV/c2 and  = (26.6 ± 2.0 ± 2.1) MeV, with √ √ c e-mail: entem@usal.es measurements at s = 4.26 GeV and s = 4.23 GeV d e-mail: fdz@usal.es 123 78 Page 2 of 10 Eur. Phys. J. C (2019) 79:78 [14]. Recently, the spin and parity of the Z c (3900) was cess one must look carefully on the analytical properties of analyzed in Ref. [15] with the result that the J P = 1+ the S-matrix. assignement is favoured by the data; moreover, its mass was In a couple of papers, Swanson [29,30] claimed that the reduced to M = (3881.2 ± 4.2 ± 52.7) and the new width Z c and Z b resonances could arise from kinematic thresh- update to  = (51.8 ± 4.6 ± 36.0) MeV, using a Flatté-like old effects, meaning that the D D̄ ∗ interaction, in case of expresion. This new determination supports the idea that the Z c (3900), is strong enough to cause an enhancement just Z c (3900)/Z c (3885) signals correspond to an unique state. above the S-wave threshold. The fact that the masses of Soon after all these experimental activities, the BESIII such resonances are slightly above their respective thresh- Collaboration reported the discovery of another charged olds seems to favor this possibility. Swanson used a simpli- state, the Z c (4020) resonance, in the e+ e− → π + π − h c fied model to describe Belle and BESIII data as kinematic channel with a mass of M = (4022.9 ± 0.8 ± 2.7) MeV/c2 cusps, successfully capturing the features of all the data and and a width of  = (7.9 ± 2.7 ± 2.6) MeV [16]. Later on, a indicating that there is no evidence for strong D D̄ ∗ or D ∗ D̄ ∗ similar structure with a mass of M = (4026.3 ± 2.6 ± 3.7) rescattering in the (I = 1) J PC = 1+− channel. Hence, MeV/c2 and a width of  = (24.8 ± 5.6 ± 5.7) MeV he concluded that isovector rescattering is not sufficiently was observed by the BESIII Collaboration in the process attractive to generate dynamical bound states and so exotic √ e+ e− → π ± (D ∗ D̄ ∗ )∓ at s = 4.26 GeV [17]. The neutral resonances are not required to explain the data. partner of the Z c (4020) was reported by the BESIII Collab- Kinematic cusp interpretations were subsequently chal- oration in Ref. [18]. lenged by Guo et al. [32]. Using a non-relativistic effective The fact that most of these structures have been found field theory, the authors argued that a cusp always exists due √ at s = 4.26 GeV led several authors to think that the to the opening of an S-wave threshold but, in order to repro- Y (4260) resonance may be responsible for the production duce the narrow and pronounced peak observed by experi- of the Z c (3900) and Z c (4020) structures due to their uncon- mentalists, a non-perturbative interaction amongst the heavy ventional properties. However, the detection of the Z c (3900) mesons is necessary and thus a nearby pole in the S-matrix at the peak of the ψ(4160) state by the CLEO-c Collaboration must appear. This motivated a detailed analysis of the exper- discarded any effect of the Y (4260) resonance [11]. imental data using reaction theory and considering different The Z c (3900) and Z c (4020) states are charged and close scenarios for the inner structure of the Z c (3900) [33]. The to D D̄ ∗ and D ∗ D̄ ∗ thresholds, respectively. Therefore, they authors claimed that a cusp is not strong enough to repro- cannot be described as pure q q̄ states and it is expected that duce the peak while a molecular or virtual state can describe D (∗) D̄ ∗ components are significant in their wave functions. equally well the experimental data. This fact can be implemented in different theoretical scenar- Further calculations with different effective field theories ios, ranging from hadron molecules [19–22], to tetraquark have been performed. Aceti et al. [34], using the local hidden structures [23–28] or simple kinematic effects linked to the gauge approach, studied the D D̄ ∗ interaction in a coupled opening of meson–meson thresholds [29,30]. channel Bethe–Salpeter calculation. They found in the I = 1 At this point, it is worth noticing that not all enhancements channel a barely D D̄ ∗ bound state decaying into ρηc and of the cross section correspond to a resonance. As discussed π J/ψ channels with a mass of 3869–3875 MeV. Addition- by Newton [31], the vicinity of the threshold of a new open- ally, Albaladejo et al. [35] performed a simultaneous study ing channel causes the cross section to show an anomalous of the invariant mass distribution for the π J/ψ and D D̄ ∗ behavior in form of a cusp. From a mathematical point of channels and found different interpretations when exploring view, this behavior is produced by the branch-point of the with an energy dependent or independent D D̄ ∗ S-wave inter- S-matrix at threshold. Resonances are also enhancements of action. In the first case, the Z c enhancement was originated the cross section, but corresponding to poles of the S-matrix from a resonance with a mass around the D D̄ ∗ threshold. In on the unphysical Riemann sheet close to the physical one. the second one, the Z c signal was produced by a virtual state In this case, the cross section can be parametrized with a which must have a hadronic molecular nature. In any case, Breit-Wigner form with two parameters: the energy (mass) a D D̄ ∗ bound state was not favored. It is worth remarking and the width. Furthermore, besides bound states which cor- that, following these authors, the tetraquark scenario would respond to poles of the S-matrix on the real energy axis of be discarded if the S-matrix pole is at the unphysical sheet the physical Riemann sheet below the lowest threshold, the (virtual state), as it would originate from D D̄ ∗ interaction S-matrix may present poles on the real energy axis but in the solely. A similar conclusion can be found in Refs. [21,36] unphysical Riemann sheet. If these virtual states have suffi- where the authors look for the Z c (3900) and Z c (4020) state ciently low energy, then their presence modifies appreciably in the π J/ψ and D D̄ ∗ and the D ∗ D̄ ∗ invariant mass spec- the energy dependence of the scattering cross section which tra respectively. The calculation was done in a quasi-potential can be anomalously large. Therefore, to correctly interpret Bethe-Salpeter equation approach at the hadronic level. They the point-wise behavior of the line shapes of a particular pro- found that the π J/ψ-D D̄ ∗ interaction produces a virtual 123 Eur. Phys. J. C (2019) 79:78 Page 3 of 10 78 state at a mass around 3870 MeV whereas they found a bound 2 Theoretical formalism state for the Z c (4020) with reasonable values of the param- eters. Additionally, the authors of Ref. [37] study the D ∗ D̄ ∗ 2.1 Constituent quark model dynamics in the local hidden gauge approach, focusing on the two pion exchange, and found that such Z c states are in The Lagrangian of Quantum Chromodynamics (QCD) with the border between weakly-bound and virtual states. massless light quarks is invariant under chiral rotations. This Many papers have been released in the last years report- symmetry does not appear in Nature indicating that it is spon- ing exploratory LQCD studies of the Z c structures [38–43]. taneously broken in QCD. Among other consequences, a Prelovsek et al. [38,39] looked for the Z c (3900)+ state in the constituent mass which depends on the quark momentum, energy region below 4 GeV, finding only two-particle π J/ψ M = M(q 2 ) and M(q 2 = ∞) = m q , is developed and and D D̄ ∗ scattering states and no signal of the Z c (3900)+ . In Goldstone-boson exchange interactions appear between the a similar search Chen et al. [40] found that the D D̄ ∗ interac- light quarks. tion is weakly repulsive and, therefore, the results do not sup- Our constituent quark model (CQM) tries to mimic port the possibility of a shallow bound state at least at the pion the previous phenomena based on the following effective mass values studied. A different approach is taken by Ikeda et Lagrangian at low-energy [65] al. [43]. They performed a coupled-channels calculation tak- ing into account the π J/ψ, ρηc and D D ∗ channels and find- L = ψ̄(i ∂/ − M(q 2 )U γ5 ) ψ, (1) ing that the interactions between them are dominated by off- diagonal couplings which may indicate that the Z c (3900)+ being U γ5 = eiλa φ γ5 / fπ the matrix of Goldstone-boson a can be explained as a threshold cusp. There are technicalities fields that can be expanded as involving LQCD computations, such as large pion masses, small volumes and a set of interpolators not large enough for i 5 a a 1 a a U γ5 = 1 + γ λ π − π π + ··· (2) having overlap with the physical state, which still prevent fπ 2 f π2 to make a definitive statement. Nonetheless, seems that the available LQCD calculations are robust enough to discard The constituent quark mass is obtained from the first term, the bound state option for the Z c+ states. the second one describes the pseudoscalar meson exchange From a molecular point of view, the fact that the Z c struc- interaction among quarks and the main contribution of the tures have I = 1 makes that the D (∗) D̄ (∗) interaction weak third term comes from the two-pion exchange which is mod- because in these channels the one-pion-interaction isospin eled by means of a scalar-meson exchange potential. coefficient is a factor 3 smaller than in the I = 0 sector. This Another nonperturbative effect is the confining interaction leaves the coupled channel calculations as the most promis- which is implemented phenomenologically in order to avoid ing option to produce some resonance, bound or virtual state, colored hadrons. In our CQM, the confinement is represented if any. as a linear potential, due to multi-gluon exchanges between In this work we analyze the I G (J PC ) = 1+ (1+− ) sector quarks, that is screened at large inter-quark distances, as a in a coupled-channels scheme, including the closest meson– consequence of sea quarks [66]: meson thresholds. The meson–meson interaction is described   in the framework of a constituent quark model1 successfully r ) = −ac (1 − e−μc r ) + VCON (  qc · λ (λ  q̄c ). (3) employed to explain the meson and baryon phenomenology from the light to the heavy quark sector [46–55]. Moreover, Here, ac and μc are model parameters. One can see that the D (∗) D̄ ∗ interaction deduced from the model has been the potential is linear at short inter-quark distances with an satisfactorily used to describe meson–meson [56–58] and effective confinement strength σ = −ac μc (λ c · λ  c ), while i j meson–baryon [59–62] molecular states such as the X (3872) it becomes constant at large distances. as a D D̄ ∗ molecule coupled to cc̄(n 3 P1 ) states [63,64]. Beyond the nonperturbative energy scale one expects that The structure of the present manuscript is organized in the the dynamics of the bound-state system is governed by QCD following way. In Sect. 2 the theoretical framework is briefly perturbative effects. We take it into account through the one- presented and discussed. Section 3 is devoted to the analysis gluon exchange potential derived from the following vertex and discussion on the obtained results. We summarize and Lagrangian give some conclusions in Sect. 4.  Lqqg = i 4π αs ψ̄γμ G aμ λa ψ, (4) where αs is the strong coupling constant, λa are the SU (3) μ 1 The interested reader is referred to Refs. [44,45] for detailed reviews colour matrices and G a is the gluon field. To consistently about the quark model in which this work is based. treat the light- and heavy-quark sectors we employ a gluon 123 78 Page 4 of 10 Eur. Phys. J. C (2019) 79:78 coupling constant that scales with the reduced mass of the where E is the energy of the system. The direct potential interacting quarks. Its explicit expression can be found in,  , P) RGM V αα ( P  of a reaction AB → C D can be written as D e.g., Ref. [46]. A detailed physical background of the quark model can be  RGM VDαα ( P , P)  = d p A d p B d pC d p D found in Refs. [46,67]. The model parameters and explicit i, j expressions for the potentials can be also found therein. We want to highlight here that the interaction terms between × φC∗ ( pC )φ ∗D ( p D )Viαα   j ( P , P)φ A ( p  A )φ B ( p B ). (7) light–light, light–heavy and heavy–heavy quarks are not the same in our formalism, i.e. while Goldstone-boson where {i, j} runs over the constituents of the involved exchanges are considered when the two quarks are light, mesons. The quark rearrangement potential RGM V Rαα ( P , P)  they do not appear in the other two configurations: light– represents a natural way to connect meson–meson channels heavy and heavy–heavy; however, the one-gluon exchange with different quark content, such as π J/ψ and D D̄ ∗ , and it and confining potentials are flavor-blind. The reason is that is given by the Goldstone-boson exchanges originates in the sponta-  neous breakdown of the Chiral Symmetry. This symmetry is RGM V αα ( P R  , P)  = d p A d p B d pC d p D d P φ ∗A ( pC ) explicitly broken at the level of heavy quarks and therefore i, j no exchanges appear associated with these quarks. However, × φ ∗D ( p D )Viαα    A )φ B ( p B )δ (3) ( P − P ) , j ( P , P )Pmn φ A ( p the presence of a light quark in the charm mesons allows to naturally incorporate the pion exchange interaction in the (8) D (∗) D̄ ∗ dynamics. where Pmn is the operator that exchanges quarks between 2.2 Resonating group method and Lippmann–Schwinger clusters. equation The meson eigenstates φC ( pC ) are calculated by means of the two-body Schrödinger equation, using the Gaussian The aforementioned CQM specifies the microscopic inter- Expansion Method [69]. This method provides enough accu- action among constituent quarks. In order to describe the racy and simplifies the subsequent evaluation of the needed interaction at the meson level we employ the Resonating matrix elements. With the aim of optimizing the Gaussian Group Method [68], where mesons are considered as quark– ranges employing a reduced number of free parameters, we antiquark clusters and an effective cluster–cluster interaction use Gaussian trial functions whose ranges are given by a emerges from the underlying q q̄ dynamics. geometrical progression [69]. This choice produces a dense We assume that the wave function of a system composed distribution at short distances enabling a better description of two mesons A and B with distinguishable quarks can be of the dynamics mediated by short range potentials. written as2 The solution of the coupled-channels RGM equations is performed deriving from Eq. (6) a set of coupled Lippmann–   p A p B P Pc.m. |ψ = φ A ( p A )φ B ( p B )χα ( P), (5) Schwinger equations of the form  where φC ( pC ) is the wave function of a general meson C with Tαα (E; p , p) = Vαα ( p , p) + dp p 2 Vαα ( p , p ) pC the relative momentum between the quark and antiquark α of the meson C. The wave function which takes into account 1 the relative motion of the two mesons is χα ( P),  where α × T α (E; p , p), (9) E − Eα ( p ) α labels the set of quantum numbers needed to uniquely define a certain partial wave. where Vαα ( p , p) is the projected potential that contains The projected Schrödinger equation for the relative wave the direct and rearrangement potentials, and Eα ( p ) is the function can be written as follows: energy corresponding to a momentum p , written in the non-   P 2  relativistic case as: − E χα ( P ) + RGM αα   VD ( P , P) 2μ p2 α Eα ( p) = + Mα . (10) 2μα + RGM V Rαα ( P , P)  χα ( P)  d P = 0, (6) Here, μα is the reduced mass of the AB system cor- 2 responding to the channel α, and Mα is the difference Note that, for simplicity of the discussion presented herein, we have dropped off the spin–isospin wave function, the product of the two color between the threshold of the AB system and the one we singlets and the wave function that describes the center-of-mass motion. take as a reference. 123 Eur. Phys. J. C (2019) 79:78 Page 5 of 10 78 We solve the coupled-channels Lippmann–Schwinger π π equation using the matrix-inversion method proposed in Ref. [70], generalized in order to include channels with different A A A thresholds. Once the T -matrix is calculated, we determine T B B the on-shell part which is directly related to the scattering B matrix (in the case of nonrelativistic kinematics): (a) (b)  Sαα = 1 − 2πi μα μα kα kα Tαα (E + i0+ ; kα , kα ), (11) Fig. 1 Background (a) and rescattering (b) contributions for the pro- duction of the two-meson AB channel with kα the on-shell momentum for channel α. √ Our aim is to explore the existence of states above and rest frame of π + A + B at energy s. The on-shell momenta below thresholds within the same formalism. Thus, we have are given by to analytically continue all the potentials and kernels for com- √ plex momenta in order to find the poles of the T -matrix in λ1/2 ( s, m AB , m π ) kπ Z c = √ , (15) any possible Riemann sheet. 2 s λ1/2 (m AB , m A , m B ) k AB = , (16) 2.3 Production line shapes of the Z c 2m AB Most of the experimental data of the Z c (3900) and Z c (4020) where λ(M, m 1 , m 2 ) = [(M 2 − m 2+ )(M 2 − m 2− )], with √ m ± = m 1 ± m 2 . We will assume that, in our phenomeno- were taken at s = 4.26 GeV, which corresponds to the energy of the resonance Y (4260), though there are several logical vertex, the production is mainly produced through √ the S-wave, which is expected to dominate near threshold as measurements at different energies such as s = 4.17 GeV [11], which may question the usual claim that the the D-wave is suppressed due to its momentum dependence. Y (4260) is the parent of the Z c (3900). To simplify the Integrating the angles, we have theoretical description of the Z c generating process and to avoid possible complexities of considering the internal 1 k AB kπ Z c d = |Mβ (m AB )|2 dm AB , (17) structure of the Y (4260) resonance or other possible ori- (2π )3 4s gin of the Z c (3900), we adopt a phenomenological vertex which creates the π + A + B triplet, where A + B is the with β the quantum numbers of the channel AB. final two-body state of the coupled-channels calculation, i.e. The Lorentz-invariant production amplitude, M, describes A + B = {π J/ψ, ρηc , D D̄ ∗ , D ∗ D̄ ∗ }. the Z c → AB production. It is diagramatically shown in The three-body decay π + A + B of a resonance in its rest Fig. 1 and can be written as frame can be written as Mβ (m AB ) ⎛ ⎞ (2π )4   β β (E + i0+ ; p, k β ) d = |M|2 d(P; pπ , p A , p B ), (12) = ⎝A β − Aβ d3 p T ⎠. 2M p 2 /2μ − E − i0+ β where M is the production amplitude and d is the three- (18) body phase-space given by where β ( ) denotes a given AB channel in the coupled-   channels calculation. We only consider production through  d 3 pi d(P; p1 , p2 , p3 ) = δ (4) P − pi . the AB in S-wave (as the vertex does not have momentum (2π )3 2E i i i dependence) and we add the amplitudes Aβ to take into (13) account different production probabilities for each channel AB. Therefore, we can express d in terms of the invariant mass Experimentalists measure events in the σ (e+ e− → spectrum of the two-meson channel, m AB , as π ± Z c∓ ) × B(Z c∓ → AB) process and thus, to describe the data, we need to add a normalization factor to translate the 1 k AB kπ Z c decay rate into events: d = |M|2 dm AB d AB dπ Z c , (14) (2π )5 16 s d Z c →AB N (m AB ) = N AB × . (19) where (k AB ,  AB ) is the on-shell momentum of the AB pair dm AB and π Z c is the solid angle of the π in the center-of-mass 123 78 Page 6 of 10 Eur. Phys. J. C (2019) 79:78 In order to fit {A AB , N AB } for each channel, we minimize 40 a global χ 2 function using all the available experimental data for channels π J/ψ, D D̄ ∗ and D ∗ D̄ ∗ : 30   N the (xi ) − N exp (xi ) 2 χ 2 ({A, N }) = exp , (20) i σi 20 and the uncertainty of the parameters {A, N } are estimated 10 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 0 3900 3950 4000 4050 4100  2 N − N BF χ (N ) ≈ χ (N ) + 2 2 BF , (21) σN 200 with σN the error of the best-fit value N BF . 150 3 Results 100 As already mentioned above, we perform a coupled-channels calculation within the framework of the constituent quark 50 model of Ref. [46] for the I G (J PC ) = 1+ (1+− ) sector, including the closest thresholds to the experimental masses of 0 the Z c (3900) and Z c (4020) hadrons, that is: π J/ψ (3234.19 3900 3950 4000 4050 4100 MeV/c2 ), ρηc (3755.79 MeV/c2 ), D D̄ ∗ (3875.85 MeV/c2 ), D ∗ D̄ ∗ (4017.24 MeV/c2 ), where the threshold masses are shown in parenthesis. The h c π channel is not taken into Fig. 2 √Line shapes for D D̄ ∗ (upper pannel) and π J/ψ (lower pan- nel) at s = 4.26 GeV. Experimental data are from Refs. [14,15], account because it couples weakly to the D (∗) D̄ ∗ channels. In respectively. The theoretical line shapes have been convoluted with the fact, the only contribution of this channel would come from experimental resolution. The line-shape’s 68% uncertainty is shown as the 3 D1 component of the internal wave function of the D ∗ a shadowed area meson, which in our model is ∼ 0.03% and it is neglected in the present calculation. For the 3 S1 component of the D ∗ , 60 the coupling to h c π is exactly zero. Similarly, we do not 50 include other nearby channels such as the χc J ρ ones, whose couplings are found to be almost three orders of magnitude 40 smaller than those of the J/ψπ and ηc ρ channels. Our results for the invariant mass distribution of the D D̄ ∗ , 30 π J/ψ and D ∗ D̄ ∗ channels are shown in, respectively, the 20 upper and lower panels of Figs. 2 and 3. Note that all chan- nels mentioned in the above paragraph are included in the 10 calculation of the theoretical line shapes. To translate the 0 decay rates into events, we use a normalization factor N AB 4020 4040 4060 4080 4100 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 Fig. 3 Line shapes for D ∗ D̄ ∗ at s = 4.26 GeV. Experimental data are from Ref. [17]. The theoretical line shapes have been convoluted the statistical 68%-confident level (CL) of the fit, obtained by with the experimental resolution. The line-shape’s 68%-uncertainty is propagating the errors of the fitted parameters by means of shown as a shadowed area the covariance matrix. The χ 2 fit is performed with the exper- imental results of Refs. [9,14,15], that is, experimental data for D D̄ ∗ and π J/ψ channels. The amplitudes A AB obtained factors and amplitudes are shown in Table 1, the result on the from such fit are, then, employed for the D ∗ D̄ ∗ channel in χ 2 /d.o.f. is also collected therein. In order to describe the order to obtain the global normalization N D ∗ D̄ ∗ from exper- experimental measurement, the theoretical line shapes have imental data of Ref. [17]. The values for the normalization been convoluted with the detector resolution. 123 Eur. Phys. J. C (2019) 79:78 Page 7 of 10 78 Table 1 Production amplitudes and normalization factors for the line shapes which are fitted using Eq. (20). The minimum value of the 200 χ 2 /d.o.f. is also given. The 68%-uncertainties in the amplitudes and normalization factors are obtained from the fit 150 Channel N AB (×107 ) A AB π J/ψ 3.76 ± 0.09 0.34 ± 0.01 100 D D̄ ∗ 0.80 ± 0.04 0.76 ± 0.01 D ∗ D̄ ∗ 19.33 ± 0.7 0.66 ± 0.01 50 ρηc (1S) 0.55 ± 0.10 −1.00 ± 0.04 0 2 /d.o.f. χmin 1.89 3900 3950 4000 4050 4100 √ Fig. 4 Line shapes for ηc ρ channel at s = 4.26 GeV (red line). D D̄ ∗ data was taken from the m For the fit, the D D̄ ∗ thresh- Experimental data, shown as blue dots, is taken from Ref. [71]. The old up to 3.92 GeV. This is the region where the signal of the theoretical line shape has been convoluted with the experimental reso- Z c (3900) dominates over the background [13], from there lution. The shadowed area shows the 68% uncertainty of the line shape 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 1.2 D D̄ ∗ channel one can clearly see two enhancements related with the opening of the D D̄ ∗ and D ∗ D̄ ∗ thresholds, which 1.0 can be associated with the Z c (3900) and the Z c (4020). In 0.8 the π J/ψ case, only one enhancement appear around 3.87 GeV while the opening of the D ∗ D̄ ∗ channel appears as a 0.6 slight step down in the number of events. Besides the unavoidable normalization factors and ampli- 0.4 tudes related to the phenomenological π + A + B vertex, our 0.2 formalism is able to reproduce the experimental data with- out fine-tuning the parameters of the constituent quark model 0.0 3900 3950 4000 4050 4100 that, in fact, are those which determine the pole position. The meson–meson interaction driven by the quark rearrangement process, which connect the D (∗) D̄ ∗ channels with the π J/ψ Fig. 5 Line shapes for different coupled-channels calculations: Only and ρηc ones, is of the same order of magnitude than the direct D D̄ ∗ channel (red, long-dashed line), D D̄ ∗ + D ∗ D̄ ∗ (blue, dashed line), interaction, D (∗) D̄ ∗ − D (∗) D̄ ∗ . This points to an important ρηc + D D̄ ∗ (green, dot-dashed line), ρηc + D D̄ ∗ + D ∗ D̄ ∗ (orange, dotted line) and π J/ψ +ρηc + D D̄ ∗ + D ∗ D̄ ∗ (black, solid line). All the role of the π J/ψ(ρηc ) − D ∗ D̄ (∗) mixings in this sector, in line shapes are normalized to coincide at the Z c (3900)’s peak position, agreement with the recent claims of some lattice studies [43], according to the normalization of Table 1. See text for discussion which concluded that the Z c cannot be described as a sim- ple D D̄ ∗ molecule. Recently, data on the ρηc channel has been published on Ref. [71]. Although this dataset has not ble with the experimental situation. On the other hand, adding been included in the global fit of A AB and N AB parame- to the D D̄ ∗ channel the D ∗ D̄ ∗ one generates a second struc- ters of Table 1, done using solely the experimental data of ture at its threshold opening associated to the Z c (4020) peak. π J/ψ and D D̄ ∗ , we present in Fig. 4 the invariant mass dis- One can see that the second enhancement is much higher than tribution of the ρηc channel. The Nρηc global normalization the experimental data if we only include D (∗) D̄ ∗ channels. is, however, adjusted from the available experimental data The line shape moves closer to the experimental situation employing those same A AB amplitudes used for all the line when the ρηc and π J/ψ channels are considered in the cal- shapes. culation, showing that they play an important role in building To disentangle the contribution of the different channels to the observed enhancements. the line shape, we compare in Fig. 5 the line shape for the full To deepen into the nature of the Z c (3900) and the calculation with those taking into account partial combina- Z c (4020), we have examined the analytic structure on the tions of the considered channels. When the D D̄ ∗ scattering complex energy plane of the S-matrix for the different channel is considered alone, a small enhancement basically coupled-channels calculations. Our results are shown in due to the pion tensor interaction between the S and D waves Table 2. For the Z c (3900), one can see that even for a one- appears but the generated peak is too wide. The inclusion of channel D D̄ ∗ calculation the S-matrix shows a pole below the ρηc channel narrows that peak, making it more compati- threshold in the second Riemann sheet. When the ρηc is 123 78 Page 8 of 10 Eur. Phys. J. C (2019) 79:78 Table 2 The S-matrix pole positions, in MeV/c2 , for different coupled- The RS columns indicate whether the pole has been found in the first channels calculations. The included channels for each case are shown (F) or second (S) Riemann sheet of a given channel. Each channel in in the first column. Poles are given in the second and fourth columns the coupled-channels calculation is represented as an array’s element, by the value of the complex energy in a specific Riemann sheet (RS). ordered with increasing energy Calculation Z c (3900) pole RS Z c (4020) pole RS D D̄ ∗ 3871.37 − 2.17 i (S) – – D D̄ ∗ + D ∗ D̄ ∗ 3872.27 − 1.85 i (S,F) 4014.16 − 0.10 i (S,S) ρηc + D D̄ ∗ 3871.32 − 0.00 i (S,S) – – ρηc + D D̄ ∗ + D ∗ D̄ ∗ 3872.07 − 0.00 i (S,S,F) 4013.10 − 0.00 i (S,S,S) π J/ψ + ρηc + D D̄ ∗ + D ∗ D̄ ∗ 3871.74 − 0.00 i (S,S,S,F) 4013.21 − 0.00 i (S,S,S,S) Table 3 Comparison of our S-matrix pole position, in MeV/c2 , with positions found in this work with other estimations from stud- several other calculations from literature ies of the Z c (3900)/Z c (4020) resonances. The pole position Calculation Z c (3900) Type described by other studies spans over a wide range of values in the complex plane, which could be caused by respective This work 3871.74 Virtual model details. However, as we can appreciate, most studies Ref. [34] 3878 − 23 i Resonance favour the virtual assignation, as in this work. Ref. [35] 3894 ± 6 ± 1 − 30 ± 12 ± 6 i Resonance 3886 ± 4 ± 1 − 22 ± 6 ± 4 i Resonance 3831 ± 26+7 −28 Virtual 4 Summary 3844 ± 19+12 −21 Virtual Ref. [43] 3709 ± 94 − 183(46) i Virtual The Z c (3900) and Z c (4020) are very peculiar structures 3748 ± 76 − 157(32) i Virtual different from other molecular resonances of the charmo- 3686 ± 56 − 44(27) i Virtual nium spectrum. As they have I = 1, the diagonal interaction Ref. [72] 3876 − 5 i Resonance between the D (∗) D̄ (∗) channels is too suppressed to develop Calculation Z c (4020) Type resonances, being the non-diagonal rearrangement interac- tion due to the coupling with other channels responsible for This work 4013.21 Virtual the structures appeared in the line shapes. Therefore, they Ref. [37] (3990–4000) − 50 i Bound/Virtual should only appear in calculations involving coupled chan- nels. To show up this mechanism, we have performed, within included a pole corresponding with a virtual state appears. the framework of a constituent quark model, a coupled- The inclusion of the rest of the channels does not change dras- channels calculation of the I G (J PC ) = 1+ (1+− ) sector tically the pole position. However, it is necessary to reproduce around the energies of the Z c (3900) and Z c (4020), includ- the experimental line shape. In the complete calculation, the ing the most relevant thresholds. The line shapes of the D D̄ ∗ , Z c (3900) is associated with a pole located in the imaginary π J/ψ and D ∗ D̄ ∗ invariant mass distributions are well repro- axis and the second Riemann sheet below the D D̄ ∗ threshold duced without any fine-tuning of the model parameters. The and, thus, it is a virtual state. The situation is similar in the analysis of the S-matrix poles allows us to conclude that the case of the Z c (4020), which is interpreted as a virtual state point-wise behavior of the line shapes is due to the presence located below the D ∗ D̄ ∗ threshold. of two virtual states that can be seen as D (∗) D̄ ∗ threshold It is worth noticing that the real part of the Z c (3900) cusps, and have an overall good agreement with the Z c (3900) mass pole, shown in Table 2, is always below the D D̄ ∗ and Z c (4020) signals. These results confirm the conclusion threshold and, hence, it does not seem to be compatible with of the lattice QCD calculation of Ref. [43]. the mass and width estimations from experimental measure- Acknowledgements This work has been partially funded by Minis- ments [9,14,15]. Obviously, in our calculation, the Z c (3900) terio de Economía, Industria y Competitividad under Contracts no. and Z c (4020) are not resonances but virtual state and, there- FPA2016-77177-C2-2-P, FPA2014-55613-P, FPA2017-86989-P and fore, a directly comparison of the complex energy of the SEV-2016-0588 and by Junta de Castilla y León and ERDF under Con- pole with the Breit–Wigner parametrization of a resonance tract No. SA041U16. P.G.O. acknowledges the financial support from Spanish MINECO’s Juan de la Cierva-Incorporación programme, Grant would be misleading. The only way to connect our results Agreement no. IJCI-2016-28525. J.S. acknowledges the financial sup- with the experimental data is through the description of the port from the European Union’s Horizon 2020 research and innovation line shapes. In Table 3 we show a comparison of the pole programme under the Marie Skłodowska-Curie Grant Agreement No. 123 Eur. Phys. J. C (2019) 79:78 Page 9 of 10 78 665919, and from Spanish MINECO’s Juan de la Cierva-Incorporación 24. J.M. Dias, F.S. Navarra, M. Nielsen, C.M. Zanetti, Phys. Rev. D programme, Grant Agreement no. IJCI-2016-30028. 88, 016004 (2013). arXiv:1304.6433 [hep- ph] 25. S.S. Agaev, K. Azizi, H. Sundu, Phys. Rev. D 96, 034026 (2017). Data Availability Statement This manuscript has no associated data arXiv:1706.01216 [hep-ph] or the data will not be deposited. [Authors’ comment: All associated 26. Z.-G. Wang, T. Huang, Phys. Rev. D 89, 054019 (2014). data have already been included in the paper.] arXiv:1310.2422 [hep-ph] 27. C.-F. Qiao, L. Tang, Eur. Phys. J. C 74, 3122 (2014). Open Access This article is distributed under the terms of the Creative arXiv:1307.6654 [hep-ph] Commons Attribution 4.0 International License (http://creativecomm 28. C. Deng, J. Ping, F. Wang, Phys. Rev. D 90, 054009 (2014). ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, arXiv:1402.0777 [hep-ph] and reproduction in any medium, provided you give appropriate credit 29. E.S. Swanson, Phys. Rev. D 91, 034009 (2015). arXiv:1409.3291 to the original author(s) and the source, provide a link to the Creative [hep-ph] Commons license, and indicate if changes were made. 30. E.S. Swanson, Int. J. Mod. Phys. E 25, 1642010 (2016). Funded by SCOAP3 . arXiv:1504.07952 [hep-ph] 31. R.G. Newton, Scattering Theory of Waves and Particles (Springer, Berlin, Heidelberg, 1982). https://doi.org/10.1007/ 978-3-642-88128-2 32. F.-K. Guo, C. Hanhart, Q. Wang, Q. Zhao, Phys. Rev. D 91, 051504 References (2015). arXiv:1411.5584 [hep-ph] 33. A. Pilloni, C. Fernandez-Ramirez, A. Jackura, V. Mathieu, M. 1. S.K. Choi et al. (Belle), Phys. Rev. Lett. 91, 262001 (2003). Mikhasenko, J. Nys, A.P. Szczepaniak (JPAC), Phys. Lett. B 772, arXiv:hep-ex/0309032 [hep-ex] 200 (2017). arXiv:1612.06490 [hep-ph] 2. B. Aubert et al. (BaBar), Phys. Rev. D 71, 071103 (2005). 34. F. Aceti, M. Bayar, E. Oset, A. Martinez Torres, K.P. Khemchan- arXiv:hep-ex/0406022 [hep-ex] dani, J.M. Dias, F.S. Navarra, M. Nielsen, Phys. Rev. D 90, 016003 3. D. Acosta et al. (CDF), Phys. Rev. Lett. 93, 072001 (2004). (2014). arXiv:1401.8216 [hep-ph] arXiv:hep-ex/0312021 [hep-ex] 35. M. Albaladejo, F.-K. Guo, C. Hidalgo-Duque, J. Nieves, Phys. Lett. 4. V.M. Abazov et al. (D0), Phys. Rev. Lett. 93, 162002 (2004). B 755, 337 (2016). arXiv:1512.03638 [hep-ph] arXiv:hep-ex/0405004 [hep-ex] 36. J. He, D.-Y. Chen, Eur. Phys. J. C 78, 94 (2018). arXiv:1712.05653 5. B. Aubert et al. (BaBar), Phys. Rev. Lett. 90, 242001 (2003). [hep-ph] arXiv:hep-ex/0304021 [hep-ex] 37. F. Aceti, M. Bayar, J.M. Dias, E. Oset, Eur. Phys. J. A 50, 103 6. D. Besson et al. (CLEO), Phys. Rev. D 68, 032002 (2003). [Erra- (2014). arXiv:1401.2076 [hep-ph] tum: Phys. Rev.D75,119908(2007)]. arXiv:hep-ex/0305100 [hep- 38. S. Prelovsek, C.B. Lang, L. Leskovec, D. Mohler, Phys. Rev. D 91, ex] 014504 (2015). arXiv:1405.7623 [hep-lat] 7. A. Bondar et al. (Belle), Phys. Rev. Lett. 108, 122001 (2012). 39. S. Prelovsek, L. Leskovec, Phys. Lett. B 727, 172 (2013). arXiv:1110.2251 [hep-ex] arXiv:1308.2097 [hep-lat] 8. I. Adachi et al. (Belle) (2012). arXiv:1209.6450 [hep-ex] 40. Y. Chen et al., Phys. Rev. D 89, 094506 (2014). arXiv:1403.1318 9. M. Ablikim et al. (BESIII), Phys. Rev. Lett. 110, 252001 (2013). [hep-lat] arXiv:1303.5949 [hep-ex] 41. Y. Chen et al. (CLQCD), Phys. Rev. D 92, 054507 (2015). 10. Z.Q. Liu et al. (Belle), Phys. Rev. Lett. 110, 252002 (2013). arXiv:1503.02371 [hep-lat] arXiv:1304.0121 [hep-ex] 42. S.-H. Lee, C. DeTar, H. Na, D. Mohler (Fermilab Lattice, MILC), 11. T. Xiao, S. Dobbs, A. Tomaradze, K.K. Seth, Phys. Lett. B 727, (2014). arXiv:1411.1389 [hep-lat] 366 (2013). arXiv:1304.3036 [hep-ex] 43. Y. Ikeda, S. Aoki, T. Doi, S. Gongyo, T. Hatsuda, T. Inoue, T. Iritani, 12. M. Ablikim et al. (BESIII), Phys. Rev. Lett. 115, 112003 (2015). N. Ishii, K. Murano, K. Sasaki (HAL QCD), Phys. Rev. Lett. 117, arXiv:1506.06018 [hep-ex] 242001 (2016). arXiv:1602.03465 [hep-lat] 13. M. Ablikim et al. (BESIII), Phys. Rev. Lett. 112, 022001 (2014). 44. A. Valcarce, H. Garcilazo, F. Fernandez, P. Gonzalez, Rep. Prog. arXiv:1310.1163 [hep-ex] Phys. 68, 965 (2005). arXiv:hep-ph/0502173 [hep-ph] 14. M. Ablikim et al. (BESIII), Phys. Rev. D 92, 092006 (2015). 45. J. Segovia, D.R. Entem, F. Fernandez, E. Hernandez, Int. J. Mod. arXiv:1509.01398 [hep-ex] Phys. E 22, 1330026 (2013). arXiv:1309.6926 [hep-ph] 15. M. Ablikim et al. (BESIII), Phys. Rev. Lett. 119, 072001 (2017). 46. J. Vijande, F. Fernandez, A. Valcarce, J. Phys. G 31, 481 (2005). arXiv:1706.04100 [hep-ex] arXiv:hep-ph/0411299 [hep-ph] 16. M. Ablikim et al. (BESIII), Phys. Rev. Lett. 111, 242001 (2013). 47. A. Valcarce, H. Garcilazo, J. Vijande, Phys. Rev. C 72, 025206 arXiv:1309.1896 [hep-ex] (2005). arXiv:hep-ph/0507297 [hep-ph] 17. M. Ablikim et al. (BESIII), Phys. Rev. Lett. 112, 132001 (2014). 48. D.R. Entem, F. Fernandez, Phys. Rev. C 73, 045214 (2006) arXiv:1308.2760 [hep-ex] 49. A. Valcarce, H. Garcilazo, J. Vijande, Eur. Phys. J. A 37, 217 18. M. Ablikim et al. (BESIII), Phys. Rev. Lett. 113, 212002 (2014). (2008). arXiv:0807.2973 [hep-ph] arXiv:1409.6577 [hep-ex] 50. J. Segovia, D.R. Entem, F. Fernandez, Phys. Lett. B 662, 33 (2008) 19. F.-K. Guo, C. Hidalgo-Duque, J. Nieves, M.P. Valderrama, Phys. 51. J. Segovia, A.M. Yasser, D.R. Entem, F. Fernandez, Phys. Rev. D Rev. D 88, 054007 (2013). arXiv:1303.6608 [hep-ph] 80, 054017 (2009) 20. T. Mehen, J. Powell, Phys. Rev. D 88, 034017 (2013). 52. P.G. Ortega, D.R. Entem, F. Fernandez, Phys. Lett. B 696, 352 arXiv:1306.5459 [hep-ph] (2011) 21. J. He, X. Liu, Z.-F. Sun, S.-L. Zhu, Eur. Phys. J. C 73, 2635 (2013). 53. J. Segovia, D.R. Entem, F. Fernandez, Proceedings, 11th Interna- arXiv:1308.2999 [hep-ph] tional Workshop on Meson Production, Properties and Interaction 22. X.-H. Liu, L. Ma, L.-P. Sun, X. Liu, S.-L. Zhu, Phys. Rev. D 90, (MESON 2010): Cracow, Poland, June 10–15, 2010. Phys. Rev. D 074020 (2014). arXiv:1407.3684 [hep-ph] 83, 114018 (2011) 23. A. Esposito, A.L. Guerrieri, F. Piccinini, A. Pilloni, A.D. Polosa, 54. J. Segovia, D.R. Entem, F. Fernandez, Phys. Rev. D 91, 094020 Int. J. Mod. Phys. A 30, 1530002 (2015). arXiv:1411.5997 [hep-ph] (2015). arXiv:1502.03827 [hep-ph] 123 78 Page 10 of 10 Eur. Phys. J. C (2019) 79:78 55. J. Segovia, P.G. Ortega, D.R. Entem, F. Fernndez, Phys. Rev. D 93, 64. P.G. Ortega, D.R. Entem, F. Fernandez, J. Phys. G 40, 065107 074027 (2016). arXiv:1601.05093 [hep-ph] (2013). arXiv:1205.1699 [hep-ph] 56. P.G. Ortega, J. Segovia, D.R. Entem, F. Fernandez, Phys. Rev. D 65. D. Diakonov, Prog. Part. Nucl. Phys. 51, 173 (2003). 94, 074037 (2016). arXiv:1603.07000 [hep-ph] arXiv:hep-ph/0212026 [hep-ph] 57. P.G. Ortega, J. Segovia, D.R. Entem, F. Fernndez, Phys. Rev. D 95, 66. G.S. Bali, H. Neff, T. Duessel, T. Lippert, K. Schilling (SESAM), 034010 (2017). arXiv:1612.04826 [hep-ph] Phys. Rev. D 71, 114513 (2005). arXiv:hep-lat/0505012 [hep-lat] 58. P.G. Ortega, J. Segovia, D.R. Entem, F. Fernndez, Phys. Lett. B 67. J. Segovia, A.M. Yasser, D.R. Entem, F. Fernandez, Phys. Rev. D 778, 1 (2018). arXiv:1706.02639 [hep-ph] 78, 114033 (2008) 59. P.G. Ortega, D.R. Entem, F. Fernandez, Phys. Lett. B 718, 1381 68. J.A. Wheeler, Phys. Rev. 52, 1083 (1937) (2013). arXiv:1210.2633 [hep-ph] 69. E. Hiyama, Y. Kino, M. Kamimura, Prog. Part. Nucl. Phys. 51, 223 60. P.G. Ortega, D.R. Entem, F. Fernndez, Phys. Rev. D 90, 114013 (2003) (2014) 70. R. Machleidt, in Computational Nuclear Physics. 2: Nuclear Reac- 61. P.G. Ortega, D.R. Entem, F. Fernndez, Phys. Lett. B 729, 24 (2014) tions, ed. by K. Langanke, J.A. Maruhn, SE koonin (Springer, 62. P.G. Ortega, D.R. Entem, F. Fernndez, Phys. Lett. B 764, 207 Berlin, 1993), pp. 1-29 (2017). arXiv:1606.06148 [hep-ph] 71. C.-Z. Yuan, Int. J. Mod. Phys. A 33, 1830018 (2018). 63. P.G. Ortega, J. Segovia, D.R. Entem, F. Fernandez, Phys. Rev. D arXiv:1808.01570 [hep-ex] 81, 054023 (2010). arXiv:0907.3997 [hep-ph] 72. J. He, Phys. Rev. D 92, 034004 (2015). arXiv:1505.05379 [hep-ph] 123