A combined analysis of the $Z_c(3900)$ and the $Z_{cs}(3985)$ exotic states

We have performed a combined analysis of the BESIII data for both the $Z_c(3900)$ and $Z_{cs}(3985)$ structures, assuming that the latter is an SU(3) flavor partner of the former one. We have improved on the previous analysis of Albaladejo $et$ $al.$ [Phys. Lett. B 755, 337 (2016)] by computing the amplitude for the $D_1\bar{D}D^*$ triangle diagram considering both $D$ and $S$-wave $D_1D^*\pi$ couplings. We have also investigated effects from SU(3) light-flavor violations, which are found to be moderate and of the order of 20%. The successful reproduction of the BESIII spectra, in both the hidden-charm and hidden-charm strange sectors, strongly supports that the $Z_{cs}(3985)$ and $Z_c(3900)$ are SU(3) flavor partners placed in the same octet multiplet. The best results are obtained when an energy-dependent term in the diagonal $D^{(*)}\bar D_{(s)}^{(*)}$ interaction is included, leading to resonances (poles above the thresholds) to describe these exotic states. We have also made predictions for the isovector $Z_{c}^*$ and isodoublet $Z_{cs}^*$, $D^*\bar{D}^*$ and $D^*\bar{D}_{s}^*$ molecules, with $J^{PC}=1^{+-}$ and $J^{P}=1^{+}$, respectively. These states would be heavy-quark spin symmetry partners of the $Z_{c}$ and $Z_{cs}$. Besides the determination of the masses and widths of the $Z_c(3900)$ and $Z_{cs}(3985)$, we also predict those of the $Z_c^*$ and $Z_{cs}^*$ resonances.


I. INTRODUCTION
The discovery of the χ c1 (3872) [2], also known as X(3872), in 2003 built a landmark in the study of strong interactions and opened the gate to the abundance of the XY Z structures in the heavy quarkonium region. Many of them are difficult to be understood from the conventional quark model point of view and thus turn out to be excellent candidates for exotic states, i.e., hadrons with a content other than a quark-antiquark pair (qq) or three quarks (qqq). A large amount of experimental and theoretical studies are devoted to those XY Z states, see e.g. Refs. [3][4][5][6][7][8][9][10][11][12][13]. Among these states, the charged Z c (3900) [14,15] and Z cs (3985) [16] states are of particular interest since they are close to the DD * and D * D s /DD * s thresholds, respectively, with widths of the same order, which suggests that the latter be a candidate to be the strange partner of the former.
The structure associated to the Z c (3900) ± was simultaneously first observed in the J/ψπ ± mass spectrum of the e + e − → J/ψπ + π − reaction by the BESIII [14] and Belle [15] collaborations. In BESIII, the e + e − center of mass (c.m.) energy was fixed to M = 4.26 GeV, while in the Belle experiment, using the initial state radiation method, the energy was taken in the Y (4260) region. The peak was confirmed in an analysis of the CLEO-c data of the J/ψπ ± mass spectrum, with a slightly lower c.m. energy of the e + e − pair, at about 4.17 GeV [17]. Further experimental evidence for the Z c (3900) ± → J/ψπ ± came from the semi-inclusive decays of b-flavored hadrons in the D0 experiment [18], with the J/ψπ + π − invariant mass also constrained around the Y (4260) mass region.
In addition, the first evidence of the neutral Z c (3900) 0 decaying into J/ψπ 0 , was reported in Ref. [17] using the CLEO-c data. Later, a similar charged structure, named as Z c (3885) + , was observed in the (DD * ) + mass distribution by the BESIII Collaboration [19,20]. The angular analysis showed that the favored spin-parity quantum numbers for this peak were J P = 1 + . The similarity of the masses and widths of the Z c (3900) + and Z c (3885) + peaks suggests a common origin for both of them.
Extensive theoretical research has been done to understand the nature of the Z c (3900) structure, providing explanations for it as a compact tetraquark [21][22][23][24][25] or as a DD * resonant or virtual molecular state [1,[26][27][28][29][30][31][32][33][34]. The possibility that the observed peaks had a kinematic origin was also discussed in Refs. [35][36][37][38]. However, it was demonstrated in Ref. [39] that a pure kinematic two-body threshold cusp, without a near-threshold pole, is unlikely to produce a narrow peak in the mass distribution of the elastic channel of the pair of open heavy-flavor mesons, which has the threshold in the vicinity of the observed peak. 1 Furthermore, it was shown in Ref. [41] that the half-maximum width of the two-body threshold cusp lineshape is proportional 2 to 1/(µa 2 0 ), where µ and a 0 are the reduced mass and the S-wave scattering length of the hadron pair, respectively. Therefore, a pronounced peak would imply a large scattering length, unless the reduced mass is very big, and consequently it points out to the existence of a near-threshold pole. Thus, the pure two-body threshold cusp effect interpretation of the XY Z states, which show up as pronounced near-threshold peaks, is not favored, and unitarity (or re-summing up the s-channel loops) will necessarily lead to a near-threshold pole with such an interaction.
However, the situation becomes more complicated when there can be triangle singularities in the same near-threshold region. A triangle singularity (TS) is a logarithmic branch point [42][43][44], which could produce a peak resembling a resonance when it is close to the physical region (for a detailed review, we refer to Ref. [11]). The manifestation of the TS in the D 1 (2420)DD * loop in the J/ψπ ± distribution, relevant to the Z c (3900) structure in the reaction e + e − → J/ψπ + π − , was immediately noticed [26,45] after the Z c (3900) discovery. It was demonstrated that the triangle diagrams need to be properly taken into account to reproduce the spectrum, while they alone can not satisfactorily describe the J/ψπ ± mass distribution. In Ref. [1], the final-state interactions for the J/ψπ-DD * coupled channels are accounted for in the re-scattering T -matrix, calculated using the proper triangle diagrams. Such a scheme provides a satisfactory description of the e + e − → Y (4260) → J/ψπ + π − and Y (4260) → π + (D 0 D * − ) data simultaneously. Two different scnearios are considered to describe the DD * interaction, and in each case only one pole is found to be identified as the Z c state [1], which implies that the Z c (3900) + [14,15] and Z c (3885) + [19,20] have the same origin and are generated by the DD * interaction. On the other hand, the JPAC collaboration concluded that the data could also be described if there is no pole in the nearthreshold region [38]. 3 A difference between the analyses carried out in Refs. [1] and [38] is that a D-wave D 1 D * π vertex was considered in the former, while an S-wave one was assumed in the 1 In the lattice simulation of the HAL QCD Collaboration, performed with pion masses ranging from 410 to 700 MeV, a pole more than 100 MeV (with a large uncertainty) below the DD * threshold was found in the J/ψπ-ηcρ-DD * coupled-channel space [37,40]. However, the predicted DD * distribution is much broader than the experimental one reported by BESIII [20]. 2 This is the case when the interaction between the hadron pair is attractive, but it is not strong enough to form a bound state. In this case, the pole is a virtual state, and the maximum of the absolute value of the amplitude is exactly at the threshold, where a cusp is produced (see also Refs. [7,10]). 3 In the fit corresponding to the "no pole" scenario (denoted as "tr."), the JPAC analysis [38] includes a penalty in the merit χ 2 function to exclude a near-threshold pole. The authors found that this scenario cannot be statistically rejected. We nevertheless point out that even in this case, the amplitude would still allow for the existence of a pole, though its position, being far from threshold, was not given in the JPAC paper. latter.
Later on, new charged charmoniumlike structures Z c (4025) ± were observed near the (D * D * ) ± threshold in the π ∓ recoil mass spectrum [46] and another charged peak, Z c (4020) ± , was reported in the h c π ± invariant mass distribution [47]. While the mass of Z c (4020) is close to that of the Z c (4025), the reported width of the Z c (4020) is larger than that of the Z c (4025), although the resonance parameters of both exotic states agree within 1.5 σ [47]. The proximity of the Z c (4020)/Z c (4025) to the D * D * threshold suggests that it could be a proper candidate for the heavy quark spin symmetry (HQSS) partner of the Z c (3900), which would imply that the spinparity of this isovector resonance would be J P = 1 + [28,48,49].
The BESIII collaboration reported in Ref. [16] the first signal of a hidden-charm resonant structure with strangeness, Z cs (3985) − , observed in the K + recoil-mass spectrum of the reaction s and D 0D * s thresholds, which are located at 3975 MeV and 3977 MeV, respectively, immediately spurred the hadronic molecular interpretation of this state [34,[50][51][52][53][54][55][56][57][58], as a strange partner of the Z c (3900). Other possible interpretations for the nature of the Z cs (3985) were also suggested, see e.g. Refs. [59][60][61][62][63][64][65][66][67]. An analysis of the D * 0D s +D 0D * s invariant mass spectra was performed in Ref. [50], in which the D * D s and DD * s S-wave interactions (J P = 1 + ) were related to that of the DD * pair in the J P C = 1 +− channel using the SU(3) lightquark flavour symmetry. There, in addition to the direct point-like production of K + D * 0D s and K + D 0D * s , a triangle diagram mechanism with the loop D * s2D * s D 0 was also considered because of the presence of a nearby TS, which can mimic the peak structure and enhance the production of near threshold molecules [11,68]. The analysis was, in principle, improved in Ref. [58] by incorporating theD * s D * /D * D * s channels and extending the analysis to the whole energy range covered by the BESIII data. Two types of solutions that describe the data almost equally well are found in Ref. [58], and both are consistent with the interpretation of the Z cs (3985) sate as an SU(3) partner of the Z c (3900). In Ref. [69], the LHCb collaboration reported, from an analysis of the B + → J/ψφK + decay amplitude, two new hidden-charm strange structures Z cs (4000) + and Z cs (4220) + decaying into J/ψK + . In particular, the Z cs (4000) + was determined to have a mass of (4003±6 + 4 −14 ) MeV, a width of (131±15±26) MeV and J P = 1 + spin-parity. While its mass is close to that of the Z cs (3985), its width is much larger than that of the Z cs (3985), which has generated a debate on whether the Z cs (3985) and the Z cs (4000) correspond to the same state [50,56] or not [70,71] (see also Ref. [72]).
In Ref. [50], the SU(3) strangeness partner of the Z c (3900) was predicted from the pole position obtained for this state in the study of Ref. [1]. In that work, fits to the K + recoil-mas distributions of the process e + e − → K + (D * 0D s + D 0D * s ), for events collected at different energy points by BESIII, were also carried out (see Figs. 3 and 4 of Ref. [50]). A difference of a few tens of MeV was found between the two mass determinations, which is largely covered by the uncertainties.
In the present work, we perform a combined analysis of the BESIII data for both the Z c (3900) and Z cs (3985) structures assuming that the Z cs (3985) is an SU(3) flavor partner of the Z c (3900). Additionally, we improve on a different aspect. In the analysis of Ref. [1], the amplitude for the D 1D D * triangle diagram considered a D 1 D * π coupling in D-wave. However, it was found in Ref. [68] that the D-wave coupling accounts only for about half of the D 1 (2420) decay width, and a D 1 D * π coupling in S-wave is also required. The inclusion of an S-wave D 1 D * π vertex could have a significant impact on the triangle diagram mechanism for the Z c (3900) peak, since it provides a sizable different background contribution to the J/ψπ invariant mass distribution. Both the S-and D-wave couplings are considered in Ref. [33] in the tree-level diagrams, however there, the triangle mechanism is not included. Moreover, Gaussian regulators are employed in Refs. [1,50] to render the integrals in the Lippmann-Schwinger equation ultraviolet (UV) finite. While it is legitimate to use such a renormalization scheme in the small momentum regime, the interaction strength could be significantly enhanced (weakened) in the energy region far below (above) the relevant channel threshold. In order to bypass this issue and test the stability of the results obtained in Refs. [1,50], we will use in this work dimensional regularization as only contact interactions are involved.
The paper is organized as follows. In Sec. II, we briefly review the contact potentials, consistent with HQSS and SU(3) light-flavour symmetries, and construct the hadron T -scattering matrix and production amplitudes taking both the point-like and triangle diagram production mechanisms into account. The parameters related to the T -matrix are determined from a best fit to the J/ψπ, DD * and D sD * + DD * s mass distributions in various scenarios and the resulting poles are investigated in Sec. III. Section IV is devoted to studying possible SU(3) flavour violation effects. Finally, the main conclusions of this work and a brief outlook are collected in Sec. V. The explicit expression of the scalar three-point loop function is relegated to Appendix A.

A. Contact interactions
The two-particle D ( * ) (s)D ( * ) (s) hadron states can be expanded in the HQSS basis |s Q ⊗ j ; J I S [74], with s Q and j the total spin of the heavy-quark subsystem and the total angular momentum of the light degrees of freedom, respectively. For S-wave states, j is just the total spin of the lightquark subsystem and the parity of the state is P = +. In addition, the total angular momentum, J, of the state is obtained by coupling s Q ⊗ j , and I and S stand for the isospin 4 and strangeness of the light degrees of freedom. For simplicity, and when it cannot be misleading, we will omit the quantum numbers of the HQSS basis elements. In what respects to the light degrees of freedom, the D (s) and D * (s) form a 1 2 ⊗ 1 2 spin multiplet. Thus in a given isospin-strangeness sector, the where the subindex J denotes the total angular momentum. HQSS implies that the strong interaction is independent of the conserved heavy quark spins in the m Q → ∞ limit. Therefore, we have s Q ⊗ 1 Ĥ I s Q ⊗ 0 = 0, and we introduce the low energy constants (LECs) for the interactions between the D ( * ) (s)D ( * ) (s) channels, where the LECs C 0,1 depend on the isospin and strangeness of the light-quark subsystem. We will come back to this point below since SU(3) flavor symmetry will reduce the number of independent LECs to only four [49]. Then the contact potentials have the form They coincide with the interactions derived from the effective Lagrangian in Ref. [58] by identifying C d = (C 0 + C 1 )/2 and C f = (C 1 − C 0 )/2, and with those obtained in Ref. [49], taking into account that in this latter case, there is a sign difference between theD * (s) states used there and those employed in this work (see Footnote 5).
The charge conjugation of the DD , D sDs , D * D * and D * sD * s S-wave meson and antimeson states is C = (−1) J . In addition, in the J = 1 sector, while the C-parity of the mesonantimeson pair D * D * can only be negative, both positive and negative C-parities can be achieved by the combinations Thus, the potentials, V J P C , for J P C = 1 +− and J P C = 1 ++ read with the basis 1 √ 2 DD * − D * D , D * D * for J P C = 1 +− and 1 √ 2 DD * + D * D for J P C = 1 ++ , respectively, which coincide with the results of Ref. [49].
So far in the discussion, we have not explicitly considered the isospin-strangeness of the light quark subsystem. As mentioned above, the LECs C d,f depend on these quantum numbers, but SU A . From the reduction 3 ⊗3 = 1 ⊕ 8, we conclude that assuming both HQSS and SU(3) flavor symmetries, there is only a total of four independent LECs, which correspond to SU(3) singlet and octet (or equivalently isoscalar and isovector) spin LECs (C 0 and C 1 or C d and C f ) [49]. This is to say, for instance, C (1) f and C (1) d . In this notation, we can construct the singlet and octet representations as: where λ i 's are the Gell-Mann matrices. It immediately follows As the SU(3) symmetry implies that the D ( * )D( * ) interaction only distinguishes the singlet and octet representations, we introduce where C (1) and C (8) rely on the heavy-quark structure of the system. Then it is straightforward to obtain (without summing over either A or B in the following expressions and for A = B) The octet representation contains an isospin triplet (S = 0, I = 1), one isospin singlet (S = 0, I = 0) and two isospin doublets (S = ±1, I = 1 2 ). It follows that in the basis D ( * )D( * ) , 8; S I M I , and hence we find that the interactions in the channels of the Z c (3900) and of its strange partner [Z cs (3985)] are identical in the SU(3) limit [50]. As a result, one finds in the heavy quark limit that for the J P = 1 + sector with the bases in the order 1 , respectively. One can make contact with the results of Ref. [49] using C and C 1b the LECs defined therein. Note that, assuming the existence of the Z c (3900) and of the Z cs (3985), then the structure of the potentials of Eq. (11) supports the existence of their D * D * and D * D * s HQSS partners, where the first one could be identified with the Z c (4025), and the second one was first predicted in Ref. [50].
The inelastic transitions between the D (s)D * (s) and J/ψπ (J/ψK, J/ψK) are also required to describe the J/ψπ mass distribution, as well as to account for the contribution from J/ψπ and J/ψK to the widths of the Z c and Z cs , respectively. 6 The S-wave J/ψπ (J/ψK, J/ψK ) system can also be expressed in the heavy-light basis J/ψπ(K, K) J=1 = |s Q = 1 ⊗ j = 0 , with the isospin and strangeness of the state determined by the Golsdtone boson (pion, antikaon or kaon) of the state. Following the same strategy discussed above, it is straightforward to find that J/ψπ (J/ψK) only couples to 1 and D * D * (s) in J = 1 and thus using the same elements of the basis as in Eq. (11), we find where we have introduced the LEC V = − 1 ⊗ 0 Ĥ I 1 ⊗ 0 / √ 2 and assumed SU(3) flavor symmetry. The direct transitions J/ψπ → J/ψπ, J/ψK → J/ψK and J/ψK → J/ψK are neglected due to the OZI suppression. The contact interaction potentials can also be constructed from an effective Lagrangian; see, e.g., Refs. [58,76].
In this work, we will focus on the Z c (3900) and Z cs (3985) resonances, which are close to the DD * and DD * s /D * D s thresholds, respectively. To reduce the number of free parameters, we neglect the D * (s)D * (s) higher coupled-channel, and we only include explicitly J/ψπ and DD * − D * D / √ 2, or 6 In addition to the D (s)D * (s) and J/ψπ(K) channels, the ηcρ(K * ) threshold is located in the energy range of interest. However, regarding the Z c(s) structures and the D (s)D * (s) and J/ψπ invariant mass distributions, neglecting the ηcρ(K * ) channel is a sounded approximation. On one hand, the ηcρ(K * ) interaction is Okubo-Zweig-Iizuka (OZI) suppressed and no significant structure is observed near its threshold. On the other hand, the ηcρ(K * ) threshold is relatively far from the Z c(s) peak and the J/ψπ(K) threshold, which makes it possible to absorb, in an effective way, its contribution in the inelastic D (s)D * (s) → J/ψπ(K) transition contact term. In principle, the channel ηcρ(K * ) could also be included explicitly, once the data of interest for this channel were statistically significant. However, currently the signal of the Zc(3900) in the ρηc channel [75] is much less significant/important than that in the J/ψπ channel. Nevertheless, it is worth exploring whether the HALQCD observation that the ηcρ channel has a strong coupling with the DD * channel [37] is consistent with this experimental fact. J/ψK and DD * s − D * D s / √ 2, labelled as channels 1 and 2, respectively. As commented above, the first channel acts as a support or decay channel. Then the couple-channel contact interactions in the Z c (3900) and Z cs (3985) sectors read where the superscript s stands for the strange sector. In what follows, we parameterize the matrix elements of the potential as 7 HQSS and SU (3) light flavor symmetries imply thatṼ s 12 =Ṽ 12 andṼ s 22 =Ṽ 22 . 8 The coupled-channel T -matrix can be obtained by where with m i,1 and m i,2 the masses of the two mesons in the ith channel. Note that the use of physical masses in the above loop functions breaks SU(3) symmetry. The loop functions are logarithmically divergent and need to be regularized. In Refs. [1,50], Gaussian form factors are introduced into the potentials to render the G i (s) function UV well defined, which however distorts the interaction strength for energies far from thresholds. In this work, we evaluate the loop function G i (s) with a once-subtracted dispersion relation and its explicit expression reads for s ≥ (m i,1 + m i,2 ) 2 [77] Re G i (s) = 1 16π 2 a i (µ) + log and µ the renormalization scale. The G i (s) function should be µ-independent, and a change of µ should be compensated by that of a i (µ). The value of a i (µ) could be estimated by matching the above expression for G i (s) to the loop function regularized by a hard cutoff q max of the order of 1-1.5 GeV (see for instance, Eqs. (51) and (52) of Ref. [78]). Taking into account that we have to evaluate the function G i (s) not only for real s ≥ (m i,1 +m i,2 ) 2 , but also below threshold and in the second Riemann sheet (RS) of the complex-s plane as well, to look for the position of resonances, we refer to Appendix A of Ref. [79] for its 7 Here the √ 2mπ factor is introduced only for dimensional arguments, and does not imply that we treat nonrelativistically the pion. 8 Note that the use of mD or mD s in Eq. (14), from the normalization pre-factors of the heavy-fields, introduces a small SU(3) breaking correction.

B. Production amplitudes
The T -matrix introduced in Eq. (15) accounts for the final-state re-scattering, and one additionally needs to construct a model for the production amplitudes for the e + e − → J/ψππ, e + e − → DD * π, and e + e − → K + D * 0D s /K + D 0D * s reactions. In Refs. [1,26,38,45,68,80], it was demonstrated that the D 1 (2420)DD * triangle diagram is important for the understanding of the Z c (3900) . It can produce a resonance-like structure, which enhances the production of near-threshold resonances [11]. In addition to the triangle-diagram production mechanism, we will also include the direct point-like production in this work. The Feynman diagrams for the e + e − → J/ψππ and DD * π are shown in Fig. 1, where the reactions proceed through the formation of the resonance Y (4230)/Y (4260). The value of the coupling Y D 1 D is irrelevant for the description of the line shapes as it can be effectively absorbed into the overall normalization factors. Regarding the D 1 D * π vertex, in Refs. [1,26] a D-wave coupling is used supported by the small width of D 1 (2420), (31.3 ± 1.9) MeV [81], which suggests that it is approximately a charmed meson with j P = 3 2 + , with j the total angular momentum of its light degrees of freedom, and thus decays into the D ( * ) π mainly in D-wave. The D-wave D 1 D * π coupling is described by the Lagrangian which respects HQSS [68,82], where F π = 92.1 MeV and σ i are the pion decay constant and the spin Pauli matrices, respectively. On the other hand, the super-fields H a and T i a , represent the j = 1 2 − and 3 2 + spin multiplets, respectively. In addition, Tr[·] denotes the trace in the spinor space, and φ ba collects the pion fields with a, b the light flavor indices The D-wave coupling h D is determined from the central value of the D 2 width to be |h D | = 1.17 GeV −1 [68], which in turn leads to 15.2 MeV for the D 1 (2420) width, that is, only about half of the total. By requiring that the decay of the D 1 is saturated by the D * π (and sequential Dππ) mode, an S-wave D 1 D * π coupling is required to account for the rest of the D 1 width [68], with |h S | = 0.57, where the pion-energy factor is introduced because the pions are pseudo-Goldstone bosons of the spontaneous breaking of chiral symmetry. We denote by A 1 (s, t) and A 2 (s, t) the amplitudes of the Y → J/ψπ + π − and Y → π + D * − D 0 decays, respectively. Here s (t) stands for the square of the invariant masses of the J/ψπ − or D * − D 0 (J/ψπ + or D * − π + ) pairs, respectively, for each of the two reactions. Then, up to some irrelevant constants, the amplitude A 1 (s, t) reads where q ± and E π ± are the three-momentum and the energy of the π ± , respectively, and the α term accounts for the D-wave contribution of the point-like production (panel (1c) in Fig. 1). The S-wave π + π − re-scattering is not explicitly included in the amplitude and its contribution will be modeled by a symmetric smooth background. On the other hand, T 12 (s) represents the scattering amplitude of DD * → J/ψπ − obtained in Eq. (15) and I(s) is the scalar three-point (D 1D D * ) loop function where P 2 = M 2 , k 2 = m 2 π , s = (P − k) 2 with M the total c.m. energy of e + e − system. The expression for the triangular loop function of Eq. (23) is explicitly given in Appendix A. After the appropriate sum and average over the polarizations, one obtains 9 where q 2 π (s) = λ(M 2 , s, m 2 π )/(4M 2 ) with λ(x, y, z) = x 2 + y 2 + z 2 − 2xy − 2yz − 2xz the Källén function, E π (s) = (M 2 + m 2 π − s)/(2M ), and θ is the relative angle between the two pions in the Y (4230)/Y (4260) rest frame. In addition, Likewise, one finds where β accounts for the π + D * − D 0 point-like production, i.e. panel (2c) of Fig. 1. In Ref. [1], only the D-wave D 1 → D * π transition was considered. To assess the impact of the S-wave vertex, we will consider in this work two cases: only D-wave, i.e. h S = 0, as in Ref. [1], and both S-and D-wave couplings.
The data on e + e − → K + (D 0D * s + D * 0D s ) are measured within the e + e − c.m. energy region M ∈ [4.268, 4.698] GeV [16]. In this energy region, especially for M = 4.681 GeV, it is found that the D * s2 (2573)D * s D 0 triangle diagrams, see Fig. 2, can facilitate the production of the Z cs (3985) resonance. The D * s2 DK vertex can be described by the SU(3) extension of the Lagrangian of Eq. (18). Since only the D-wave coupling is involved for D * s2 → D ( * ) K, its strength is irrelevant as it can be reabsorbed into the normalization factor for line shapes. For the e + e − energy values measured by BESIII, it is natural to assume that the reaction e + e − → K + (D 0D * s + D * 0D s ) proceeds through the ψ(4660) resonance, e + e − → ψ(4660) → K + (D 0D * s + D * 0D s ). We denote by where q 2 K (s) = λ(M 2 , s, m 2 K )/(4M 2 ), G P 1 P 2 is the loop function of Eq. (17) evaluated with P 1 and P 2 the particles running in the two-body loop, and I s (s) represents the D * s2D * s D scalar triangle integral. Here we have used the relation and r represents the point-like production of K + D 0D * s /K + D * 0D s and is a parameter accounting 10 The contributions of the production through D * s2 and the point-like one in Ref. [50] should be added incoherently, instead of coherently as done in Ref. [50], because their accompanied K + are in D-and S-waves, respectively. The impact on the numerical results of Ref. [50] is marginal due to the very small contribution from the interference terms.
For the Y → J/ψπ + π − and Y → π + D * − D 0 reactions, the J/ψπ − and D 0 D * − spectra are obtained from the amplitudes as where N i is an unknown normalization factor, t ± i (s) are the limits of the Mandelstam variable t for the decay mode i (see, e.g., the Dalitz plot section of the Kinematics review in Ref. [81]), and B i (s) is an incoherent background mimicking possible contributions from crossed channels, misidentified events, and higher waves other than the S-wave. For the J/ψπ − and D 0 D * − spectra, the background is parameterized as [1,14] with m 1,− = m J/ψ + m π , m 2,− = m D * + m D , and m + = M − m π . In this work, we will fix d i = 1 to reduce the number of free parameters. We have checked that releasing d i as a free parameter barely improves the fit quality.
For the e + e − → K + (D 0D * s + D * 0D s ), we will be only interested in the energy region close to the D * D s and DD * s thresholds and we will only fit to the low-energy tail. Hence, we will simply use the background (B) employed in the experimental analysis of Ref. [16]. The line shape for the K + recoil mass distribution at the different total e + e − energy values M can be computed as [50] with t ± s,i (s) the limits of the Mandelstam variable t for the decay mode i. The integrated luminosity L int , the detection efficiency¯ and the correction factor f corr can be found in Ref. [16], and m ψ = 4630 MeV and Γ ψ = 62 MeV are employed for the mass and width of the ψ(4660) resonance, respectively [81]. Finally, the factor N s is associated with the e + e − annihilation vertex, and here we take the same value for the e + e − energy range from 4.628 to 4.698 GeV.

III. DESCRIPTION OF THE EXPERIMENTAL SPECTRA: THE SU(3) FLAVOUR SYMMETRIC CASE
In this section, we will determine the hidden-charm two-meson scattering T -matrix of Eq. (15), assuming SU(3) flavour symmetry, from a combined fit to the J/ψπ − , D 0 D * − , and the K + recoil-mass [RM(K + )] distributions of the e + e − → J/ψπ + π − at c.m. energies of 4.23 and 4.26 GeV [73], e + e − → π + D 0 D * − at 4.26 GeV [20], and e + e − → K + (D * 0 D − s + D 0 D * − s ) at 4.628, 4.641, 4.661, 4.681 and 4.698 GeV [16]. The leading-order (LO) approximation for the contact potentials of Eq. (13) amounts to approximating them simply by constants. However, such interaction kernels cannot generate resonances above the higher threshold, even when coupled channel effects are considered. To avoid this problem and to take into account the possibility of the Z c(s) poles being located above the corresponding open-charm thresholds, we introduce an energy-dependent term inṼ (s) 22 , and fixṼ (s) 12 to a constant, since the J/ψπ (J/ψK) threshold is located far from the Z c (Z cs ),Ṽ (s) where C Z , b, and C 12 are constants. The LO constant potential case corresponds to b = 0, and we will consider both scenarios below, i.e., b = 0 and b = 0. Before proceeding to describe the data, we need to discuss the subtraction constants a i (µ) in the G i (s) functions defined in Eq. (17), which can be determined by matching the loop function evaluated at threshold with a hard cutoff or a Gaussian form factor with the cutoff Λ around 1 GeV [78,83]. We have taken a i (µ) = a s i (µ), consistent with SU(3) flavor symmetry. Nevertheless, it is easy to check that variations of a 1 (µ) can be fully absorbed into changes of the C Z , b, and C 12 LECs. However, this is not the case for the subtraction constant a 2 (µ). 11 Thus, in this work we have fixed a 1 (µ) = −2.77, with µ = 1 GeV, to match the G(s) function at the J/ψπ threshold to that in Ref. [1], regulated using a Gaussian form factor with Λ = 1.5 GeV. For a 2 (µ = 1 GeV), we have considered two different values, −2.5 and −3.0, although only the results for a 2 (µ) = −3.0 will be shown in plots. The LECs obtained from fits with a 2 (µ) = −2.5 will be included in the uncertainties of the pole analysis and are collected in Table II. As commented above, an alternative way to evaluate the loop integral in Eq. (16) is to use a hard cutoff q max , see e.g. Ref. [85]. The subtraction constant a 2 (µ) may be estimated by comparing with the loop function obtained using a hard cutoff (q max ) regularization, with a natural value for q max . Matching the loop function using the two regularizations evaluated at the DD * threshold, a 2 (µ = 1 GeV) = −2.5 and −3.0 correspond to q max = 1.2 and 1.8 GeV, respectively.
In summary, we have only three 12 (C Z , C 12 , b) free parameters for the T -matrix, four normalization factors (N 1 and N 1 for e + e − → J/ψπ + π − at 4.26 and 4.23 GeV, respectively, N 2 for e + e − → π + D 0 D * − at 4.26 GeV, and a common one N s for the e + e − → ψ(4660) → K + (D * 0 D − s + D 0 D * − s ) reaction for all e + e − c.m. energies), three background parameters (see Eq. (30)) for each of the e + e − → Y (4230) → J/ψπ + π − , e + e − → Y (4260) → J/ψπ + π − and e + e − → Y (4260) → π + D * − D 0 , and three additional parameters α, β and r related to the point-like production of J/ψπ + π − , π + D 0D * and K + (DD ( * ) s + D * D s ), respectively. There is a total of 19 real (or 18 in the case of constant potentials) fit parameters, to be determined from fitting to more than 220 data points. Adding a relative phase between the TS and the point-like production mechanisms barely improves the fits.
In this work, we will consider four different fit schemes, which correspond to considering or neglecting the strength (h S ) of the S-wave D 1 D * π vertex and the energy-dependent term of the V (s) 22 diagonal part of the interaction between the two heavy-light mesons (parameter b in Eq. (32)), • Scheme IA: We use only the D-wave D 1 D * π coupling (h S = 0) and take constant potentials (b = 0).
• Scheme IB: We extend the Scheme IA by considering also the energy-dependent part of the two heavy-light meson interaction (b = 0).
• Scheme IIB: We extend the Scheme IIA by considering also the energy-dependent part of the two heavy-light meson interaction (b = 0).
As mentioned above, for each scheme two different values for the subtraction constant, a 2 (µ = 1 GeV) = −2.5 and −3.0 are considered. Results of the four fit scenarios are collected in Table I, where only the parameters of the two-meson T -matrix are compiled. The results obtained for constant potentials (b = 0), i.e., in Schemes IA and IIA, are shown in Figs. 3 and 4, for the Z c and Z cs structures, respectively. The two fits lead to similar χ 2 /dof, and the best fit curves are barely different, in particular, for the J/ψπ invariant mass distributions. However, they correspond to distinctly different backgrounds [cf. Eq. (30)], especially for the Y (4230) → J/ψππ reaction. The difference can be traced back to that the D-wave D 1 D * π coupling would induce D-wave pions, which lead to a different Dalitz plot projection onto the J/ψπ invariant mass distribution compared with the S-wave ones; see Ref. [68]. Nonetheless, the scattering T -matrix is well constrained.
As found in Ref. [1] and shown here in Table I, the fit quality is significantly improved once the strength (b) of the energy-dependent potential term is allowed to vary (Schemes IB and IIB). The spectra found for these two fits are shown now in Figs. 5 and 6 for the hidden-charm and hidden-charm strange sectors, respectively. Both fits provide similarly quite good descriptions of the experimental distributions, with only tiny differences. As in the case of Schemes IA and IIA, displayed Figs. 3 and 4, the inclusion of the S-wave D 1 D * π coupling modifies the relative weight of the background in the non-strange hidden charm reactions.
A full discussion on the spectroscopic content of our fits will be given below, but we can anticipate here our main conclusion. We see that the BESIII data sets related to the Z c (3900)/Z c (3885) and Z cs (3985) signatures can be simultaneously reproduced using SU(3) flavour symmetry, which supports the assumption that the Z cs (3985) is the strange partner of the Z c (3900)/Z c (3885). Yet, the Z cs structure is only significant in the data set at M = 4.681 GeV, likely due to the enhancement induced by the triangle singularity from the diagram (3b) in Fig. 2, as already discussed in Refs. [50,58]. Data with higher statistics will be highly valuable.
Resonances and virtual states are identified as poles in different RSs of the T -matrix. For a twochannel problem, one has a 4-sheet Riemann surface in the complex energy plane. The four RSs can be accessed through the analytical continuation of the G i (s) loop functions. The expression given in Eq. (17) stands for the physical or first RS, namely, G i (s) = G I i (s). The loop function in the unphysical RS, G II i (s), which is continuously connected to G I i (s) through the cut spanning from the threshold of the ith channel to infinity along the positive real-s axis, is obtained by analytic continuation as [79]  , which is the sign of the imaginary part of the c.m. momentum of particles in the ith channel, the four RSs can be labeled as RS ±± . In this convention, the physical RS for the 2-coupled channels is labeled as RS-I = RS ++ , and the other three RSs are: RS-II = RS −+ , RS-III = RS −− , and RS-IV = RS +− . RS-II is connected to the physical region through the interval between the 1st and 2nd thresholds, and RS-III is connected to the physical region above the 2nd threshold along the real axis. RS-IV is not directly connected to the physical region, however the poles located on RS-IV can still leave impact on the physical observables. The real and imaginary parts of a pole are identified as the mass and half width, respectively, for the corresponding resonance.
For Schemes IB and IIB, the poles are found above the D * D (s) threshold on RS-III, as shown in Table II. They correspond to resonances which can decay to J/ψ/π (J/ψK) and DD * (D * D s /DD * s ). 13 The results for the non-strange sector are in agreement within errors with those found in Ref. [1] for the Λ 2 = 1.0 GeV case. We postpone the comparison of the predictions found here for the strange sector with those obtained in our previous analysis of Ref. [50] to the next section, where some SU(3) breaking effects will be included.
For Schemes IA and IIA (b = 0), the description of the experimental spectra is not as good as for Schemes IB and IIB, with a larger χ 2 /dof. The poles are found located below the D * D (s) thresholds in RS-IV with negligible imaginary parts which are caused by the small C 12 off-diagonal couplings (see Table I). These poles would move into the real axis on the unphysical RS of the D * D (s) singlechannel scattering amplitude if the J/ψπ (J/ψK) channel was switched off. Therefore, the poles obtained in Schemes IA and IIA are identified as D * D (s) virtual states, which do not correspond to spatially localized particles in the sense that their spatial wave functions are not normalizable (neither does that of a resonance). However, such singularities, if located near threshold, can significantly modify the line shapes at the vicinity of the D * D (s) threshold. While the pole positions hardly change when the the subtraction constant a 2 (µ) varies from −2.5 to −3.0, the inclusion of the S-wave D 1 D * π vertex, in addition to the D-wave coupling, produces a bigger impact on the pole positions. It is interesting to notice that the real part of the virtual state pole in Schemes IA  Table I; so are the bands in the following plots. The energy resolution in the BESIII measurement is very high [73], much finer than the bin size, and can be neglected. Thus, the effects due to the finite 15 MeV bin size for the e + e − → J/ψππ reaction (upper panels) are accounted for by averaging over the energies of each bin.
and IIA is about 100 MeV below the corresponding threshold, in line with the HAL QCD result for the Z c (3900) [37,40]. This occurs when the attraction is relatively weak. However, when the pole is located above threshold, as in Schemes IB and IIB, it is much closer to threshold.
In the heavy quark limit, the interactions of the isovector J P C = 1 +− D * D * and of the isodoublet J P = 1 + D * D * s pairs are equal to those of DD * and DD * s − D * D s / √ 2, respectively (see Eqs. (11) and (12)). The existence of the Z c and Z cs hints at the possible existence of HQSS partner resonances with the same quantum numbers. Based on the LECs obtained from the different fits compiled in Table I  masses and half widths of the Z * c and Z * cs are collected in Table II. FIG. 6. Same as in Fig. 4, but for Schemes IB and IIB, which incorporate a non-vanishing energy-dependent part of the potential (b = 0).  The spin partners predicted in Table II were obtained by neglecting the channel coupling between the pseudoscalar-vector and vector-vector meson pairs. Yet, there is a possibility that the channel coupling might be strong. The analysis done in Ref. [58] fits to the BESIII data of e + e − → K + (D * 0 D − s + D 0 D * − s ) in the whole range of RM(K + ) considering constant contact interactions in the DD * s /D * D s -D * D * s coupled channels. Depending on how the interaction is organized in the coupled channels, radically different fits were obtained. The fits with weak DD * s /D * D s -D * D * s channel coupling (fit 1 and fit 1 therein) agree well with that in Ref. [50]. The pole closest to the physical region is a virtual state pole, and thus also consistent with the analysis here in schemes with b = 0; in this case, spin partners of the Z c (3900) and Z cs (3985) were also predicted close to the D * D * and D * D * s thresholds, respectively. However, for fit 2 in Ref. [58], the DD * s /D * D s -D * D * s channel coupling is stronger than the diagonal interaction, and it is mainly because of the channel coupling that the Z c (3900) and Z cs (3985) are generated. In this case, the Z c (3900) (Z cs (3985)) does not have a spin partner with the same J P C (J P ). A comprehensive analysis of all data for the Z c (3900), Z cs (3985) and Z c (4020)/Z c (4025) would be helpful to finally pin down the role of the channel coupling between the pseudoscalar-vector and vector-vector meson pairs.
It is observed in Table II that the mass of the Z * c is close to that of the Z c (4020/4025), suggesting that this latter (observed) resonance is the HQSS partner of the Z c (3900). Although the Z * c width is larger than that of the Z c (4020/4025) reported by BESIII, a direct comparison is not very meaningful since the BESIII analysis was made using a Breit-Wigner parametrization. Also, as mentioned in Ref. [50], just before the Z cs (4220) peak there is a dip around the D * D * s threshold, and it is well possible that the dip is a consequence of the Z * cs (for a general discussion of the appearance of a dip-like near-threshold structure, see Ref. [41]). It is plausible that the Z cs (4000) reported by LHCb [69] has the same origin as the pole part of the Z cs (3985) signal reported by BESIII despite that the reported width of the former is much larger than that of the latter. In order to reach a firm conclusion on this, it would be very valuable to conduct in the future a joint analysis of both data sets within the framework derived here. Nevertheless, the above claim, proposed in the noted added of Ref. [50], finds support in Ref. [56], where both the BESIII and LHCb data can be described well within the same model, and with the Z cs and Z * cs appearing as virtual states.
Finally, we briefly discuss the production of the Z c(s) . In Eq. (22), the point-like production of the J/ψππ (with ππ in D-wave) is modeled by a constant coupling α, and the final state interaction for J/ψπ → J/ψπ direct transition is neglected, since it is OZI suppressed. Therefore, the pointlike production can not generate a peak structure in the J/ψπ mass distributions around the DD * threshold.

IV. SU(3) FLAVOR SYMMETRY BREAKING EFFECTS
Up to now, we have employed two-meson potentials in the SU(3) limit, which allows to use common LECs (C Z , C 12 and b) for the DD * and D * D s /DD * s interactions, cf. Eq. (32). However, the SU (3) flavor is not an exact symmetry and we could expect violations of around 20%.
In addition, we should note that the parameters related to the T -matrix are mainly determined by the Z c related mass distributions, since they incorporate more data points with smaller uncertainties, which strongly constrain the fits. Although both the Z c and Z cs , can be well described in the SU(3) limit, as demonstrated in Sec. III, cf. Figs. 5 and 6, we observe that the description of the Z cs structure at M = 4.681 GeV can be probably improved by including SU(3) breaking effects. In Table II, we observe that the mass differences between the Z cs and Z c states in Schemes IB and IIB are 98-99 MeV, evaluated using the central values, which are consistent with the predictions of Ref. [50] in the case of the SU(3) limit, i.e., 100-102 MeV evaluated using the central values obtained with different UV cutoffs. To investigate effects from light-flavor symmetry violations, we consider two ways to include SU(3) breaking terms: • (b): introducing a breaking term in the diagonal DD * s − D * D s / √ 2 potential, C s Z = C Z + δ Z , and keeping the other terms of the two-body interactions unchanged.
• (c): introducing a breaking term in the off-diagonal J/ψK-DD * s − D * D s / √ 2 potential, C s 12 = C 12 + δ 12 , and keeping the other terms of the two-body interactions unchanged.
An SU(3) breaking correction in the energy-dependent term (b-term) is not considered because this LEC is of higher order in the heavy-meson momentum expansion. On the other hand, as shown in the previous section, the pole positions are insensitive to the value of the subtraction constant a 2 (µ), thus in this section we set a 2 (µ) = −3.0 for definiteness. In Table III we show the results of these fits, and the poles so obtained are shown in Table IV. The line shapes for the Z c structure are almost unchanged after introducing the SU(3) breaking terms. The description of the RM(K + ) distributions, especially for M = 4.681 GeV, is significantly improved in schemes of type (b). However, for schemes of type (c), the description of the RM(K + ) distributions has been less improved (see Table III). Therefore and to illustrate the effects, we show the comparison between the results obtained in the SU(3) limit and the type (b) SU(3) breaking scenario in Fig. 7. In this figure, we consider a full S + D-wave D 1 D * π vertex.
While the introduction of δ 12 for the strange sector does barely improve the fit quality, it, however, leads to significantly larger uncertainties for the Z ( * ) cs parameters, see those for the schemes of type (c) in Table IV. For the schemes of type (b), the decreased value of χ 2 /dof is mainly caused by the improved description of the Z cs structure at M = 4.681 GeV. It is easy to conclude from Fig. 7 that the Z cs pole position approaches to the D * D s threshold, as can be confirmed in Table IV, which produces a more pronounced enhancement of the RM(K + ) distribution in the region close to threshold. Let us focus on Schemes IB and IIB with SU(3) breaking of type (b), which lead to smaller χ 2 /dof than Schemes IA and IIA. The mass and width of Z c are almost unchanged compared to those found in the SU(3) limit, while the mass and width of the Z cs resonance are reduced by around 15 MeV and 10 MeV, respectively. For convenience, in Fig. 8 we collect all the resonance poles in Tables II and IV. The results of this work, when SU(3) breaking effects are considered, compare reasonably well with those obtained in Ref. [50] from straight fits to the BESIII RM(K + ) distributions.
With the parameters in Table III, we estimate the size of the SU(3) breaking terms as which is consistent with a naive expectation. We stress that the ∆ b value is not the unique evaluation of the level of SU(3) violation since the LECs C Z and δ Z are not physical observables and they are renormalization scale-dependent. However, this quantity still provides a rough assessment of the SU(3) violation.

V. CONCLUSIONS
We have performed a combined analysis of the BESIII data for both the Z c (3900) and Z cs (3985) structures, and it is found that the data can be well described assuming that the latter is an SU(3) flavor partner of the former one. We have improved on the previous analysis of Ref. [1] by computing the amplitude for the D 1D D * triangle diagram considering both D-and S-wave D 1 D * π couplings. Including the S-wave D 1 D * π vertex has a certain impact on the triangle diagram mechanism for the Z c (3900) peak, since it provides a sizable different background contribution to the J/ψπ invariant mass distribution. Moreover, we have used in this work dimensional regularization to render the integrals in the Lippmann-Schwinger equation UV finite. This resolves the issue of employing Gaussian regulators, as in Refs. [1,50], which could significantly enhance (weaken) the interaction  . The fitted energy region is shaded. The energy-dependent part of potential, controlled by the LEC b is zero (different to zero and fitted to data) in the left (right) panel. In all cases, the full S + D-wave D 1 D * π vertex is used. The error bands are statistical, propagated from the uncertainties quoted in Table III. strength in the energy region far below (above) the relevant channel threshold. Finally, we have investigated effects from SU(3) light-flavor violations, which are found to be moderate and of the order of 20%. The successful reproduction of the BESIII measured spectra, in both non-strange and strange hidden-charm sectors, strongly supports that the Z cs (3985) and Z c (3900) are SU(3) flavor partners placed in the same octet multiplet. The best results are obtained when an energy-dependent term in the diagonal D ( * )D ( * ) (s) interaction is included, leading to resonances (poles above the corresponding open-charm thresholds) to describe these exotic states, though the data are also compatible with the Z c (3900) and Z cs (3985) as virtual states (poles below the corresponding open-charm thresholds on the real energy axis of the unphysical RS). We have also made predictions for the isovector Z * c and isodoublet Z * cs , D * D * and D * D * s molecules, with J P C = 1 +− and J P = 1 + , respectively. These states would be HQSS partners of the Z c and Z cs . The masses and widths of the Z c (3900), Z cs (3985), Z * c , and Z * cs resonances are collected in Table IV.   Tables II and IV, obtained in different fits studied in this work, with an energy-dependent term in the diagonal D ( * )D ( * ) (s) interaction. The respective nearby thresholds are also shown.
One important feature of the contributions from triangle diagrams with a TS close to the physical region is the sensitivity to the kinematic variables [11]. In the problem under study, the significance of the Z cs (3885) signal at the e + e − c.m. 4.681 GeV relative to the other energies can be attributed to this effect. At this respect and to better understand the nature of the Z c (3900), high-statistic data in a sufficiently large range of e + e − c.m. energies, including not only 4.23 and 4.26 GeV, but also for instance 4.29 GeV, where the TS plays a more important role, and other energies far from 4.29 GeV will be highly valuable. In this way, one should have enough information to map out the relative important on the relevant invariant mass spectra of the TS and the pole contributions.