1 Introduction

After the observation of X(3872) at Belle, more and more exotic hadrons are observed in experiment [1]. How to understand the origin and internal structure of such states becomes an important central issue of hadron physics. With the accumulation of the experimental information, the molecular state becomes one of the most important pictures to explain these exotic states. For example, the hidden charm pentaquarks observed recently at LHCb are just a little below corresponding thresholds and exhibit an excellent spectrum of S-wave molecular states [2,3,4,5,6,7,8,9].

In the tetraquark sector, the first exotic state X(3872) was soon explained as a molecular state due to its closeness to the \(D{\bar{D}}^*\) threshold [10,11,12]. In the past decade, more and more structures were observed near the \(D^{(*)}{\bar{D}}^{(*)}/B^{(*)}{\bar{B}}^{(*)}\) thresholds, such as the \(Z_b(10610)\), \(Z_b(10650)\), \(Z_c(3900)\), and \(Z_c(4020)\) [13,14,15] (see Fig. 1). Recently, BESIII Collaboration announced observation of a new exotic candidate \({Z_{cs}}(3985)^-\) in \(e^+\) \(e^-\) \(\rightarrow \) \({K^+} (D_s^-D^{*0} + D_s^{*-}D^0)\) near the \({D_{s}}^{-}{D}^{*0} / {D_{s}}^{*-}{D}^{0}\) threshold. It has a mass of \(3982.5^{+1.8}_{-2.6}\pm 2.1\,\hbox {MeV}\) and a width of \(12.8^{+5.3}_{-4.4}\pm 3.0\,\hbox {MeV}\) [16]. The LHCb also reported a structure \(Z_{cs}(4000)\) at \(4003\pm 6^{+4}_{-14}\,\hbox {MeV}\), which is very close to the \(Z_{cs}(3985)\), but has a much broader width of \(131\pm 15\pm 26\,\hbox {MeV}\) [17]. In 2009, the Y(4140) state which is below the \(D^*_s{\bar{D}}^*_s\) threshold was reported by the CDF Collaboration in the exclusive \(B\rightarrow KJ/\psi \phi \) decay [18, 19], and confirmed at D0 and CMS in 2013 [20, 21].

Since the observation of these states, much theoretical efforts have been made to understand their origin and internal structures. Besides interpretation of the X(3872) as a \(D{\bar{D}}^*\) molecular state with quantum numbers \(I^G(J^{PC})=0^-(1^{++})\), two \(Z_b\) states are also widely taken as the candidates of the \(B{\bar{B}}^*\) and \(B^*{\bar{B}}^*\) molecular states [22, 23]. The \(Z_c(3900)\) and \(Z_c(4020)\) are good partners of the \(Z_b\) states, and also be interpreted as molecular states of systems \(D{\bar{D}}^*\) and \(D^*{\bar{D}}^*\) [24,25,26,27,28]. However, As shown in Fig. 1, though these states are very close to the thresholds, the two \(Z_b\) and two \(Z_c\) states, especially the \(Z_c(3900)\), are actually beyond the thresholds, which seems to conflict with the definition of the molecular state as a loosely bound state. Moreover, many lattice calculations disfavor existence of a \(D{\bar{D}}^*\) bound state [29,30,31,32]. Therefore, many authors suggested that the \(Z_c(3900)\) is a virtual state instead of a bound state from the \(D{\bar{D}}^*\) interaction [33, 34]. Since the virtual state locates on the second Reimann sheet, its effect on the invariant mass spectrum should appear beyond the threshold. As discussed in Ref. [35], the attraction between two hadrons will exhibits itself as a structure even if it is not bound. Hence, the molecular state interpretation of these states, especially if we categorize virtual state as a molecular state, is still the most reasonable and widely accepted [36].

Fig. 1
figure 1

The experimental masses of X(3872), \(Z_c(3900)\), \(Z_c(4020)\), \(Z_b(10610)\), \(Z_c(10650)\), \(Z_{cs}(3985)\), and Y(4140) with relevant thresholds. The values are from the Review of Particle Physics [1]

The observation of the hidden charm state with strangeness \({Z_{cs}}(3985)^-\) provides more confidence of the molecular interpretation of these charmonium-like states. Since the mass of the newly observed state is nearly 100 MeV larger than \({Z_{c}}(3900)\), which is the mass difference between \({D_{s}}^{*}\) and \({D}^{*}\) meson, the state is recognized as a molecular state of \({D_{s}}^{-}{D}^{*0} + {D_{s}}^{*-}{D}^{0}\) in the literature [37,38,39,40,41,42,43,44,45]. In Ref. [37], the mass and width were well reproduced with an assignment of \({Z_{cs}}(3985)\) as a U/V-spin partner of \({Z_c}(3900)\). The states \({Z}^*_{cs}\), \(Z_{bs}\) and \({Z}^*_{bs}\) were also predicted as the U/V-spin partner states of the charged \({Z_c}(4020)\), \({Z_b}(10610)\) and \({Z_b}(10650)\), respectively. Similar conclusion can be found in Ref. [38], and the study was also extended to the \(D^{(*)}_s{\bar{D}}^{(*)}_s\) systems [39]. In Ref. [40], within a \({D_{s}}^{-}{D}^{*0}\) molecular state picture, authors found a pole located at \(3982.34 -i0.53\,\hbox {MeV}\) on the third Reimann sheet of the complex energy plane by solving the Bethe–Salpeter equation with the on-shell factorization. The structure of Y(4140) was also explained as a \({D}^{*}_{s}{{\bar{D}}}^{*}_{s}\) molecular state [46,47,48,49,50,51]. In Ref. [46], the authors proposed that the Y(4140) can be taken as a partner with hidden-strangeness of the Y(3930). In Refs. [47,48,49], the scalar \({D}^{(*)}_{s}{{\bar{D}}}^{(*)}_{s}\) molecular state was studied with the QCD sum rules, and their result favors the assignment of the Y(4140) as a molecular state. Chen et al. suggested that the Y(4140) state should be considered as a mixed state of three pure molecule states \({D}^{*0}{{\bar{D}}}^{*0}\), \({D}^{*+}{{D}}^{*-}\) and \({D}^{*+}_{s}{{D}}^{*-}_{s}\) [50]. However, most molecular state interpretations suggested spin parity \(0^+\), which is inconsistent with the quantum numbers determined at LHCb, \(J^{PC}=1^{++}\) [52]. The LHCb result disfavors the assignment of the Y(4140) as an S-wave \({D}^{(*)}_{s}{{\bar{D}}}^{(*)}_{s}\) molecular state. Anyway, the previous studies strongly suggests existence of \({D}^{(*)}_{s}{{\bar{D}}}^{(*)}_{s}\) molecular states.

In our previous studies, the \(Z_c\) and \(Z_b\) states were studied within the one-boson-exchange model by solving the Schrödinger equation or the quasipotential Bethe–Salpeter equation (qBSE) [22, 27, 53]. If we extend such study to the \(Z_{cs}\) state in the molecular state picture, a problem will appear. Different from the \(D^{(*)}{\bar{D}}^{(*)}/B^{(*)}{\bar{B}}^{(*)}\) interaction, no light meson can be exchanged between a charmed meson and a charmed strange meson. If we do not introduce the heavy meson exchange, no bound states can be produced from corresponding interactions. In Ref. [27], the \(J/\psi \) exchange was included in the study of \(D^{(*)}{\bar{D}}^{(*)}\) system, and found essential to reproduce the \(Z_c(3900)\) as in other approach [26]. And a further study of the systems \(D^{(*)}{\bar{B}}^{(*)}/B^{(*)}{\bar{D}}^{(*)}\) and \(D^{(*)}{B}^{(*)}/B^{(*)}{D}^{(*)}\) suggest that the \(J/\psi \) or \(\varUpsilon \) exchange is really important to reproduce the experimentally observed states \({Z_{c}}(3900)\), \({Z_{c}}(4020)\), \({Z_{b}}(10610)\), and \({Z_{b}}(10650)\) [54]. The new experimental observation of \({Z_{cs}}(3985)\) inspires us to explore the possible molecular states from interactions \(D^{(*)}{{\bar{D}}}^{(*)}_{s}\) / \(B^{(*)}{{\bar{B}}}^{(*)}_{s}\), \({D}^{(*)}_{s}{{\bar{D}}}^{(*)}_{s}\) / \({{B}}^{(*)}_{s}{{\bar{B}}}^{(*)}_{s}\), \({D}^{(*)}{D_{s}}^{(*)}\)/\({B}^{(*)}{B_{s}}^{(*)}\) and \({D_{s}}^{(*)}{D_{s}}^{(*)}\)/\({B_{s}}^{(*)}{B_{s}}^{(*)}\) with the \(J/\psi \) or \(\varUpsilon \) exchange included. In the calculation, the potential of the interactions will be constructed with the meson exchange and inserted into the qBSE to obtain the scattering amplitudes. The molecular states are searched for as the poles of the scattering amplitudes.

This article is organized as follows. After introduction, we present the details of theoretical frame in Sect. 2, which includes flavor wave functions, effective Lagrangians, construction of potential and a brief introduction of the qBSE approach. The numerical results for the states produced from the interaction considered will be given in Sect. 3. Finally, article ends with summary and discussion in Sect. 4.

2 Theoretical frame

First, we should construct the flavor wave functions. Since the systems considered in the current work contain one u/d quark at most, the isospins of states are fixed. We first present the wave functions for systems \(D{{\bar{D}}}^*_{s}\) and \(D^*{{\bar{D}}}^*_{s}\). Since the thresholds of these two systems are very close, under the SU(3)\(_F\) symmetry, the flavor wave functions for charged and neutral states can be expressed as [37],

$$\begin{aligned} |X_{D^*{{\bar{D}}}_{s}+D{{\bar{D}}}^*_{s}}^-\rangle&=\frac{1}{\sqrt{2}} \big (|{D}^{*0}D_{s}^{-}\rangle +\eta |{D}^{0}D_{s}^{*-}\rangle \big ),\nonumber \\ |X_{D^*{{\bar{D}}}_{s}+D{{\bar{D}}}^*_{s}}^0\rangle&=\frac{1}{\sqrt{2}} \big (|{D}^{*+}D_{s}^{-}\rangle +\eta |{D}^{+}D_{s}^{*-}\rangle \big ), \end{aligned}$$
(1)

where \(\eta =\pm \) corresponds to the \(Z_{cs}\) state with \(G_{U/V}=\pm \). To make the calculation practicable, the masses of \(D_s\) and D, as well as \(D^*_s\) and \(D^*\) mesons, will be chosen as their average value. To check the effect of such treatment, we will also consider the direct wave functions without mixing as

$$\begin{aligned} |X_{{{D}}^*{\bar{D}}_s}^-\rangle&=|{D}^{*0}D_{s}^{-}\rangle ,\quad |X_{D{{\bar{D}}}^*_{s}}^-\rangle =|{D}^{0}D_{s}^{*-}\rangle ,\nonumber \\ |X_{{{D}}^*{\bar{D}}_s}^0\rangle&=|{D}^{*+}D_{s}^{-}\rangle ,\quad |X_{D{{\bar{D}}}^*_{s}}^0\rangle =|{D}^{+}D_{s}^{*-}\rangle . \end{aligned}$$
(2)

For the \({D_{s}}{\bar{D}}_{s}^{*}\) states, the flavor function with a definite charge parity C is constructed as

$$\begin{aligned} |X_{{D_{s}}{{\bar{D}}}^*_{s}}^0\rangle&=\frac{1}{\sqrt{2}}\big (|D_{s}^{-} D_{s}^{*+}\rangle -C|D_{s}^{*-}D_{s}^{+}\rangle \big ). \end{aligned}$$
(3)

Different from the SU(3)\(_F\) symmetry, the C parity is well defined because the masses of \(D_s^+\) and \(D_s^-\), as well as the masses of \(D^{*+}_s\) and \(D^{*-}_s\) mesons, are almost the same.

The flavor wave functions for hidden charm systems \(D^*{{\bar{D}}}^*_{s}\) and \({D_{s}}{{\bar{D}}}_{s}\) can be written easily as

$$\begin{aligned} |X_{D^*{{\bar{D}}}^*_{s}}^-\rangle&=|{D}^{*0}D_{s}^{*-}\rangle , \quad |X_{D_{s}{\bar{D}}_{s}}^0\rangle =|D_{s}^{-}D_{s}^{+}\rangle ,\nonumber \\ |X_{D^*{{\bar{D}}}^*_{s}}^0\rangle&=|{D}^{*+}D_{s}^{*-}\rangle . \end{aligned}$$
(4)

For doubly charm system with strangeness \(D^*{D}_{s}/DD^*_s\), the closeness of the total masses of two channels also lead to strong mixing. Here, we construct the wave functions analogous to the \(D^*{{\bar{D}}}_{s}/D{\bar{D}}^*_s\) as

$$\begin{aligned} |X_{D^*{D}_{s}+D{D}^*_{s}}^{+}\rangle&=\frac{1}{\sqrt{2}}\big (|{D}^{*0}D_{s}^{+} \rangle +\eta |{D}^{0}D_{s}^{*+}\rangle \big ),\nonumber \\ |X_{D^*{D}_{s}+D{D}^*_{s}}^{++}\rangle&=\frac{1}{\sqrt{2}}\big (|{D}^{*+}D_{s}^{+} \rangle +\eta |{D}^{+}D_{s}^{*+}\rangle \big ). \end{aligned}$$
(5)

Though this construction is not the standard U/V spin, we still define a \(G'=\pm \) which corresponds \(\eta =\pm \). Direct wave function without mixing will be also considered as hidden charmed system \(D^*{{\bar{D}}}_{s}/D{\bar{D}}^*_s\).

The wave functions for \(D^{*}{D_{s}}^{*}\) and \(D{D_s}\) states can be constructed as,

$$\begin{aligned} |X_{D^{*}D_{s}^{*}}^{+}\rangle&=|D^{*0}D_{s}^{*+}\rangle , \quad |X_{D{D_{s}}}^{+}\rangle =|D^{0}D_{s}^{+}\rangle ,\nonumber \\ |X_{D^*D_s^*}^{++}\rangle&=|D^{*+}D_{s}^{*+}\rangle , \quad |X_{D{D_{s}}}^{++}\rangle =|D^{+}D_{s}^{+}\rangle . \end{aligned}$$
(6)

For the doubly charm system with double strangeness \(D_s^{(*)}D_s^{(*)}\), the wave function has a form of

$$\begin{aligned} |X_{D_s^{(*)}D_s^{(*)}}^{++}\rangle&=|D_{s}^{(*)+}D_{s}^{(*)+}\rangle . \end{aligned}$$
(7)

The wave functions for the hidden bottom and doubly bottom systems can be obtained analogously.

To achieve the potential kernel for the interactions, the meson exchange model will be adopted. The Lagrangians under the chiral symmetry and heavy quark limit will be introduced to depict the vertices of the meson exchanges. The Lagrangians for heavy mesons interacting with light mesons read [55,56,57,58,59],

$$\begin{aligned} {\mathcal {L}}_{{\mathcal {P}}^*{\mathcal {P}}^*{\mathbb {P}}}&= -i\frac{2g}{f_\pi }\varepsilon _{\alpha \mu \nu \lambda } v^\alpha {\mathcal {P}}^{*\mu }_{b}{{\mathcal {P}}}^{*\lambda \dag }_{a} \partial ^\nu {}{\mathbb {P}}_{ba}\nonumber \\&\quad +i \frac{2g}{f_\pi }\varepsilon _{\alpha \mu \nu \lambda } v^\alpha \widetilde{{\mathcal {P}}}^{*\mu \dag }_{a}\widetilde{{\mathcal {P}}}^{*\lambda }_{b} \partial ^\nu {}{\mathbb {P}}_{ab},\nonumber \\ {\mathcal {L}}_{{\mathcal {P}}^*{\mathcal {P}}{\mathbb {P}}}&=- \frac{2g}{f_\pi }({\mathcal {P}}^{}_b{\mathcal {P}}^{*\dag }_{a\lambda }+ {\mathcal {P}}^{*}_{b\lambda }{\mathcal {P}}^{\dag }_{a})\partial ^\lambda {} {\mathbb {P}}_{ba}\nonumber \\&\quad +\frac{2g}{f_\pi }(\widetilde{{\mathcal {P}}}^{*\dag }_{a\lambda }\widetilde{{\mathcal {P}}}_b+ \widetilde{{\mathcal {P}}}^{\dag }_{a}\widetilde{{\mathcal {P}}}^{*}_{b\lambda })\partial ^\lambda {}{\mathbb {P}}_{ab}, \nonumber \\ {\mathcal {L}}_{\mathcal {PP}{\mathbb {V}}}&= -\sqrt{2}\beta {}g_V{\mathcal {P}}^{}_b{\mathcal {P}}_a^{\dag } v\cdot {\mathbb {V}}_{ba} +\sqrt{2}\beta {}g_V\widetilde{{\mathcal {P}}}^{\dag }_a \widetilde{{\mathcal {P}}}^{}_b v\cdot {\mathbb {V}}_{ab},\nonumber \\ {\mathcal {L}}_{{\mathcal {P}}^*{\mathcal {P}}{\mathbb {V}}}&=- 2\sqrt{2}\lambda {}g_V v^\lambda \varepsilon _{\lambda \mu \alpha \beta } ({\mathcal {P}}^{}_b{\mathcal {P}}^{*\mu \dag }_a + {\mathcal {P}}_b^{*\mu }{\mathcal {P}}^{\dag }_a) (\partial ^\alpha {}{\mathbb {V}}^\beta )_{ba}\nonumber \\&\quad - 2\sqrt{2}\lambda {}g_V v^\lambda \varepsilon _{\lambda \mu \alpha \beta } (\widetilde{{\mathcal {P}}}^{*\mu \dag }_a\widetilde{{\mathcal {P}}}^{}_b +\widetilde{{\mathcal {P}}}^{\dag }_a\widetilde{{\mathcal {P}}}_b^{*\mu }) (\partial ^\alpha {}{\mathbb {V}}^\beta )_{ab},\nonumber \\ {\mathcal {L}}_{{\mathcal {P}}^*{\mathcal {P}}^*{\mathbb {V}}}&= \sqrt{2}\beta {}g_V {\mathcal {P}}_b^{*}\cdot {\mathcal {P}}^{*\dag }_a v\cdot {\mathbb {V}}_{ba}\nonumber \\&\quad -i2\sqrt{2}\lambda {}g_V{\mathcal {P}}^{*\mu }_b{\mathcal {P}}^{*\nu \dag }_a (\partial _\mu {} {\mathbb {V}}_\nu - \partial _\nu {}{\mathbb {V}}_\mu )_{ba}\nonumber \\&\quad -\sqrt{2}\beta g_V \widetilde{{\mathcal {P}}}^{*\dag }_a\widetilde{{\mathcal {P}}}_b^{*} v\cdot {\mathbb {V}}_{ab}\nonumber \\&\quad -i2\sqrt{2}\lambda {}g_V\widetilde{{\mathcal {P}}}^{*\mu \dag }_a\widetilde{{\mathcal {P}}}^{*\nu }_b(\partial _\mu {} {\mathbb {V}}_\nu - \partial _\nu {}{\mathbb {V}}_\mu )_{ab}, \nonumber \\ {\mathcal {L}}_{\mathcal {PP}{\mathbb {V}}}&= -\sqrt{2}\beta {}g_V{\mathcal {P}}^{}_b{\mathcal {P}}_a^{\dag } v\cdot {\mathbb {V}}_{ba} +\sqrt{2}\beta {}g_V\widetilde{{\mathcal {P}}}^{\dag }_a \widetilde{{\mathcal {P}}}^{}_b v\cdot {\mathbb {V}}_{ab},\nonumber \\ {\mathcal {L}}_{{\mathcal {P}}^*{\mathcal {P}}{\mathbb {V}}}&=- 2\sqrt{2}\lambda {}g_V v^\lambda \varepsilon _{\lambda \mu \alpha \beta } ({\mathcal {P}}^{}_b{\mathcal {P}}^{*\mu \dag }_a + {\mathcal {P}}_b^{*\mu }{\mathcal {P}}^{\dag }_a) (\partial ^\alpha {}{\mathbb {V}}^\beta )_{ba}\nonumber \\&\quad - 2\sqrt{2}\lambda {}g_V v^\lambda \varepsilon _{\lambda \mu \alpha \beta } (\widetilde{{\mathcal {P}}}^{*\mu \dag }_a\widetilde{{\mathcal {P}}}^{}_b +\widetilde{{\mathcal {P}}}^{\dag }_a\widetilde{{\mathcal {P}}}_b^{*\mu }) (\partial ^\alpha {}{\mathbb {V}}^\beta )_{ab},\nonumber \\ {\mathcal {L}}_{{\mathcal {P}}^*{\mathcal {P}}^*{\mathbb {V}}}&= \sqrt{2}\beta {}g_V {\mathcal {P}}_b^{*}\cdot {\mathcal {P}}^{*\dag }_a v\cdot {\mathbb {V}}_{ba}\nonumber \\&\quad -i2\sqrt{2}\lambda {}g_V{\mathcal {P}}^{*\mu }_b{\mathcal {P}}^{*\nu \dag }_a (\partial _\mu {} {\mathbb {V}}_\nu - \partial _\nu {}{\mathbb {V}}_\mu )_{ba}\nonumber \\&\quad -\sqrt{2}\beta g_V \widetilde{{\mathcal {P}}}^{*\dag }_a\widetilde{{\mathcal {P}}}_b^{*} v\cdot {\mathbb {V}}_{ab}\nonumber \\&\quad -i2\sqrt{2}\lambda {}g_V\widetilde{{\mathcal {P}}}^{*\mu \dag }_a \widetilde{{\mathcal {P}}}^{*\nu }_b(\partial _\mu {} {\mathbb {V}}_\nu - \partial _\nu {}{\mathbb {V}}_\mu )_{ab},\nonumber \\ {\mathcal {L}}_{\mathcal {PP}\sigma }&= -2g_s{\mathcal {P}}^{}_b{\mathcal {P}}^{\dag }_b\sigma -2g_s\widetilde{{\mathcal {P}}}^{}_b\widetilde{{\mathcal {P}}}^{\dag }_b\sigma ,\nonumber \\ {\mathcal {L}}_{{\mathcal {P}}^*{\mathcal {P}}^*\sigma }&= 2g_s{\mathcal {P}}^{*}_b\cdot {}{\mathcal {P}}^{*\dag }_b\sigma +2g_s\widetilde{{\mathcal {P}}}^{*}_b\cdot {}\widetilde{{\mathcal {P}}}^{*\dag }_b\sigma , \end{aligned}$$
(8)

where the velocity v should be replaced by \(i\overleftrightarrow {\partial }/2\sqrt{m_im_f}\) with the \(m_{i,f}\) being the mass of the initial or final heavy meson. \({{\mathcal {P}}}^{(*)T} =(D^{(*)0},D^{(*)+},D_s^{(*)+})\) or \((B^{(*)-},{\bar{B}}^{(*)0},{\bar{B}}_s^{(*)0})\), and satisfy the normalization relations \(\langle 0|{{\mathcal {P}}}|{Q}{\bar{q}}(0^-)\rangle =\sqrt{M_{\mathcal {P}}}\) and \(\langle 0|{{{\mathcal {P}}}}^*_\mu |{Q}{\bar{q}}(1^-)\rangle = \epsilon _\mu \sqrt{M_{{\mathcal {P}}^*}}\). The \(\mathbb P\) and \({\mathbb {V}}\) denote the pseudoscalar and vector matrices

$$\begin{aligned} {{\mathbb {P}}}= & {} \left( \begin{array}{ccc} \frac{\sqrt{3}\pi ^0+\eta }{\sqrt{6}}&{}\pi ^+&{}K^+\\ \pi ^-&{}\frac{-\sqrt{3}\pi ^0+\eta }{\sqrt{6}}&{}K^0\\ K^-&{}{\bar{K}}^0&{}-\frac{2\eta }{\sqrt{6}} \end{array}\right) ,\nonumber \\ {\mathbb {V}}= & {} \left( \begin{array}{ccc} \frac{\rho ^0+\omega }{\sqrt{2}}&{}\rho ^{+}&{}K^{*+}\\ \rho ^{-}&{}\frac{-\rho ^{0}+\omega }{\sqrt{2}}&{}K^{*0}\\ K^{*-}&{}{\bar{K}}^{*0}&{}\phi \end{array}\right) . \end{aligned}$$
(9)

The parameters involved here were determined in the literature as \(g=0.59\), \(\beta =0.9\), \(\lambda =0.56\,\hbox {GeV}^{-1}\), and \(g_s=0.76\) with \(f_\pi =132\,\hbox {MeV}\) [60,61,62,63].

In Refs. [26, 27], contribution from the \(J/\psi \) exchange is found important in the \(D{\bar{D}}^*\) interaction to produce the \(Z_c(3900)\) observed at BESIII. In the current work, we also consider such exchange with the couplings of heavy-light charmed mesons to \(J/\psi \), which are written with the help of heavy quark effective theory as [59, 64],

$$\begin{aligned} {{{\mathcal {L}}}}_{D^*_{(s)}{\bar{D}}^*_{(s)}J/\psi }= & {} -ig_{D^*_{(s)}D^*_{(s)}\psi }\big [\psi \cdot {\bar{D}}^*\overleftrightarrow {\partial }\cdot D^*\nonumber \\&- \psi ^\mu {{\bar{D}}}^* \cdot \overleftrightarrow {\partial }^\mu {D}^* +\psi ^\mu {\bar{D}}^*\cdot \overleftrightarrow {\partial } D^{*\mu } ) \big ], \nonumber \\ {{{\mathcal {L}}}}_{D_{(s)}^*{\bar{D}}_{(s)}J/\psi }= & {} g_{D^*_{(s)}D_{(s)}\psi } \, \, \epsilon _{\beta \mu \alpha \tau } \partial ^\beta \psi ^\mu \nonumber \\&({\bar{D}} \overleftrightarrow {\partial }^\tau D^{* \alpha }+{\bar{D}}^{* \alpha } \overleftrightarrow {\partial }^\tau D) , \nonumber \\ {{{\mathcal {L}}}}_{D_{(s)} {\bar{D}}_{(s)}J/\psi }= & {} ig_{D_{(s)}D_{(s)}\psi } \psi \cdot {\bar{D}}\overleftrightarrow {\partial }D, \end{aligned}$$
(10)

where the couplings are related to a single parameter \(g_2\) as

$$\begin{aligned} \frac{g_{D^*D^*\psi } }{m_{D^*}}= \frac{g_{D_{(s)}D_{(s)}\psi }}{m_D}= g_{D^*_{(s)}D_{(s)}\psi }= 2 g_2 \sqrt{m_\psi }, \end{aligned}$$
(11)

with \(g_2={\sqrt{m_\psi }}/({2m_Df_\psi })\) and \(f_\psi =405\,\hbox {MeV}\). For the bottom mesons, analogous Lagrangians can be obtained under the heavy quark symmetry for \(\varUpsilon \) exchange. The parameters can be obtained by similar relation in Eq. (11) by replacing the mass by these of bottom mesons with \(f_\varUpsilon =715\,\hbox {MeV}\) [65].

With the help of standard Feynman rule, the vertices can be easily obtained from the above Lagrangians. The potential interaction can be constructed as [9],

$$\begin{aligned} {{{\mathcal {V}}}}_{{\mathbb {P}},\sigma }= & {} I^{d,c}_i\varGamma _1\varGamma _2 P_{{\mathbb {P}},\sigma }f(q^2),\nonumber \\ {{{\mathcal {V}}}}_{{\mathbb {V}}}= & {} I^{d,c}_i\varGamma _{1\mu }\varGamma _{2\nu } P^{\mu \nu }_{{\mathbb {V}}}f(q^2), \end{aligned}$$
(12)

where the propagators are defined as usual as

$$\begin{aligned} P_{{\mathbb {P}},\sigma }= \frac{i}{q^2-m_{{\mathbb {P}},\sigma }^2},\ \ P^{\mu \nu }_{\mathbb {V}}=i\frac{-g^{\mu \nu }+q^\mu q^\nu /m^2_{{\mathbb {V}}}}{q^2-m_{\mathbb {V}}^2}. \end{aligned}$$
(13)

We introduce a form factor \(f(q^2)=\varLambda _e^2/(q^2-\varLambda _e^2)\) to reflect the off-shell effect of exchanged meson with \(m_e\) being the \(m_{{\mathbb {P}},{\mathbb {V}},\sigma }\) and q being the momentum of the exchanged meson. The current form factor can avoid overestimation of the contribution of \(J/\psi \) exchange.

In our approach, we collect the coefficients for the interaction of states as flavor factors \(I^{d,c}_i\). As discussed in Refs. [53, 54], for the hidden heavy state, the cross diagram appears in the \([D{{{\bar{D}}}^{*}_{s}}]\)/\([B{{{\bar{B}}}^{*}_{s}}]\) and \({{D}}_{s}{{\bar{D}}}^*_{s}\)/\({{B}}_{s}{{\bar{B}}}^*_{s}\) cases due to the coupling between the two parts as shown in Eqs. (1) and (3). For the states \(D{{\bar{D}}_{s}}/B{{\bar{B}}_{s}}\), \(D^{*}{{{\bar{D}}}^{*}_{s}}/B^{*}{{{\bar{B}}}^{*}_{s}}\), \(D_{s}{{\bar{D}}_{s}}/B_{s}{{\bar{B}}_{s}}\) and \(D_{s}^{*}{{{\bar{D}}}^{*}_{s}} / B_{s}^{*}{{{\bar{B}}}^{*}_{s}}\), there is no cross diagram. For the doubly heavy states, there is no cross diagram from the coupling in the \({D}^{(*)}{D_{s}}^{(*)}\)/\({B}^{(*)}{B_{s}}^{(*)}\) and the \({D_{s}}^{(*)}{D_{s}}^{(*)}\)/\({B_{s}}^{(*)}{B_{s}}^{(*)}\) cases. However, in such cases, the u channel is allowed, which provides cross diagram. In Table 1, flavor factors \(I^d_i\) and \(I^c_i\) of certain meson exchange i of certain interaction are listed for direct and cross diagrams, respectively.

Table 1 The isospin factors \(I_i^d\) and \(I_i^c\) for direct and cross diagrams and different exchange mesons. Here we use notations \({{\mathcal {P}}}^{(*)} =D^{(*)}\) or \(B^{(*)}\), \({{\mathcal {P}}}^{(*)}_s =D^{(*)}_s\) or \(B^{(*)}_s\), \([{\mathcal {P}}{\bar{{\mathcal {P}}}}_{s}^{*}]={\mathcal {P}}^*{\bar{{\mathcal {P}}}}_{s}+{\mathcal {P}}{\bar{{\mathcal {P}}}}^*_{s}\), \([{\mathcal {P}}{{{\mathcal {P}}}}_{s}^{*}]={\mathcal {P}}^*{{{\mathcal {P}}}}_{s}+{\mathcal {P}}{{{\mathcal {P}}}}^*_{s}\)

In the above, we construct the potential of the interactions considered in the current work. The scattering amplitude can be obtained with the qBSE [27, 34, 53, 66, 67]. After the partial-wave decomposition, the qBSE can be reduced to a 1-dimensional equation with a spin-parity \(J^P\) as [27],

$$\begin{aligned} i{{{\mathcal {M}}}}^{J^P}_{\lambda '\lambda }(\mathrm{p}',\mathrm{p})&=i{{{\mathcal {V}}}}^{J^P}_{\lambda ',\lambda }(\mathrm{p}',\mathrm{p})+\sum _{\lambda ''}\int \frac{\mathrm{p}''^2d\mathrm{p}''}{(2\pi )^3}\nonumber \\&\quad \cdot i{{{\mathcal {V}}}}^{J^P}_{\lambda '\lambda ''}(\mathrm{p}',\mathrm{p}'') G_0(\mathrm{p}'')i{{{\mathcal {M}}}}^{J^P}_{\lambda ''\lambda }(\mathrm{p}'',\mathrm{p}), \end{aligned}$$
(14)

where the \({{{\mathcal {M}}}}^{J^P}(\mathrm{p}',\mathrm{p})\) is partial-wave scattering amplitude, and the \(G_0(\mathrm{p}'')\) is reduced propagator with the spectator approximation in the center-of-mass frame as [27],

$$\begin{aligned} G_0= & {} \frac{\delta ^+(p''^{~2}_h-m_h^{2})}{p''^{~2}_l-m_l^{2}} =\frac{\delta ^+(p''^{0}_h-E_h)}{2E_h[(W-E_h)^2-E_l^{2}]}. \end{aligned}$$
(15)

As required by the spectator approximation, the heavier particle (remarked with h) satisfies \(p''^0_h=E_{h}(\mathrm{p}'')=\sqrt{ m_{h}^{~2}+\mathrm p''^2}\). The \(p''^0_l\) for the lighter particle (remarked as l) is then \(W-E_{h}(\mathrm{p}'')\) with W being the center-of-mass energy of the system. Here and hereafter, we define the value of the momentum in center-of-mass frame as \(\mathrm{p}=|{{\varvec{p}}}|\).

The partial wave potential can be obtained from the potential in Eq. (12) as

$$\begin{aligned} {{{\mathcal {V}}}}_{\lambda '\lambda }^{J^P}(\mathrm{p}',\mathrm{p})&=2\pi \int d\cos \theta ~[d^{J}_{\lambda \lambda '}(\theta ) {{{\mathcal {V}}}}_{\lambda '\lambda }({{\varvec{p}}}',{{\varvec{p}}})\nonumber \\&\quad +\eta d^{J}_{-\lambda \lambda '}(\theta ) {{{\mathcal {V}}}}_{\lambda '-\lambda }({{\varvec{p}}}',{{\varvec{p}}})], \end{aligned}$$
(16)

where \(\eta =PP_1P_2(-1)^{J-J_1-J_2}\) with P and J being parity and spin for system, and two constituent heavy mesons. The initial and final relative momenta are chosen as \({{\varvec{p}}}=(0,0,\mathrm{p})\) and \({{\varvec{p}}}'=(\mathrm{p}'\sin \theta ,0,\mathrm{p}'\cos \theta )\). The \(d^J_{\lambda \lambda '}(\theta )\) is the Wigner d-matrix. An exponential regularization was also introduced as a form factor into the reduced propagator as \(G_0(\mathrm{p}'')\rightarrow G_0(\mathrm{p}'')e^{-2(p''^2_l-m_l^2)^2/\varLambda _r^4}\) [27].

3 Numerical results

With the preparation above, we can obtain the scattering amplitudes of the interaction of hidden heavy and doubly heavy systems, and the molecular states can be searched for as the poles in the complex energy plane. Because we are only interesting in the pole of the scattering amplitude, we only need to find the position where \(|1-V(z)G(z)|=0\) with z being the complex continuation of the center-of-mass energy of the system W [27]. In addition, we take two free parameters \(\varLambda _e\) and \(\varLambda _r\) as \(\varLambda \) for simplification.

3.1 Numerical results with single-channel calculation

In Tables 2, 3, 4 and 5, we present the results for the binding energy \(E_B\) with single-channel calculation. Here, the binding energy is defined as \(E_B=M_{th}-z\), with the \(M_{th}\) and z being the threshold and the position of the pole of the bound state. For single-channel calculation, the pole is at real axis of the complex energy plane. In the current work, we only consider the spin parities which can be produced in S-wave. Here we scan the values of cutoff \(\varLambda \) in a range smaller than 5 GeV and present results with some selected values of cutoff \(\varLambda \) if there is a bound state produced from the corresponding interaction with a binding energy smaller than 40 MeV.

The results for the systems of two pseudoscalar mesons are listed in Table 2. Only state with spin parity \(0^+\) can be produced in S-wave for two pseudoscalar mesons. In the charmed sector, the bound state can be produced from the hidden charmed systems \(D{\bar{D}}_s\) and \(D_s{\bar{D}}_s\). However the cutoff required to produce the former state is much larger than the later one. For the hidden bottom systems \(B{\bar{B}}_s\) and \(B_s{\bar{B}}_s\), the interactions are still attractive, but a cutoff larger than 5 GeV is required to produce a \(B{\bar{B}}_s\) bound state. The attraction in the \(B{\bar{B}}_s\) interaction is strong enough to form a bound state at a cutoff of about 1 GeV. No bound states is found in the doubly charmed and doubly bottom systems, \(DD_s\), \(D_sD_s\), \(BB_s\), and \(B_sB_s\). Hence, our results support the existence of deeply bound states \(D_s{\bar{D}}_s\) and \(B_s{\bar{B}}\) with \(0^+\), and the \(D{\bar{D}}_s\) interaction are considerably attractive.

Table 2 The binding energies \(E_B=M_{th}-z\) of the bound states from the interactions \(D_{(s)}{\bar{D}}_s/B_{(s)}{\bar{B}}_s\) and \(D_{(s)}{\bar{D}}_s/B_{(s)}{\bar{B}}_s\) with some selected values of cutoff \(\varLambda \). The “\(\cdot \cdot \cdot \)” means that no bound state is found in the considered range of the cutoff. The cutoff \(\varLambda \) and binding energy W are in the units of GeV and MeV

The results for the systems of a pseudoscalar meson and a vector meson are listed in Table 3. The spin parities for these systems are \(1^+\) in S wave. Here, the bound states with different \(G_{U/V}\) spins introduced in Eq. (1) are listed for the hidden charmed strange system \(D^*{{\bar{D}}}_{s}+D{\bar{D}}^*_s\). The calculation supports the existence of a state with \(G_{U/V}=+\) at a cutoff of about 3 GeV, which is relevant to the state \(Z_{cs}(3985)\) observed at BESIII recently. No bound state can be found with \(G_{U/V}=-\) from the \(D^*{{\bar{D}}}_{s}+D{\bar{D}}^*_s\) interaction. In the bottom sector, the \(B^*{{\bar{B}}}_{s}+B{\bar{B}}^*_s\) interaction with \(G_{U/V}=+\) is still attractive while a cutoff larger than 5 GeV is required to produce a bound state. There is still no bound state produced from the interaction with \(G_{U/V}=-\). For the doubly charmed system, we still consider the system with \(G'\), a \(D^*{{D}}_{s}+DD_{s}^*\) state with \(G'=-\) can be found at a cutoff of about 2 GeV, and a state with the same \(G'\) spin is produced from the interaction \(B^*{{B}}_{s}+BB_{s}^*\) at a smaller cutoff of about 0.8 GeV.

Table 3 The binding energies of the bound states from the interactions \(D_{(s)}{\bar{D}}^*_s/B_{(s)}{\bar{B}}^*_s\), and \(D_{(s)}{D}^*_s/B_{(s)}{B}^*_s\) with some selected values of cutoff \(\varLambda \). The \(G_{U/V}\), \(G'\), and C are defined in Eqs. (1, 3, 5). Other notations are the same as Table 2

The \(D^*{\bar{D}}_s+D{\bar{D}}^*_s\) bound state with \(G_{U/V}=+\) can be related to the \(Z_{cs}(3895)\) state observed at BESIII. In the literature, it is suggested as the U/V-spin partner of the \(Z_c(3900)\) observed before the \(Z_{cs}(3985)\) [37]. It should be reminded that light meson exchanges have considerable contribution to binding of \(D{\bar{D}}^*\) system. In our previous study, the calculation suggests that the \(J/\psi \) exchange is essential to reproduce the \(Z_c(3900)\) [27, 54]. As shown in Table 1, only \(J/\psi \) meson involves in \(D^{(*)}{{\bar{D}}}^{(*)}_{s}\) interaction. Hence, the \(J/\psi \) exchange is essential to reproduce both \(Z_{cs}(3985)\) and \(Z_c(3900)\). In Table 4, we present the binding energies of the \(D^{(*)}{\bar{D}}^{(*)}\) states with only \(J/\psi \) exchange and of the \(D^{(*)}{{\bar{D}}}^{(*)}_{s}\) systems. Generally speaking, there is no significant difference between two cases as expected, which means it is reasonable to consider \(D^{(*)}{{\bar{D}}}^{(*)}_{s}\) state as strange partner of \(D^{(*)}{{\bar{D}}}^{(*)}\) state in the SU(3)\(_F\) symmetry. Moreover, we can expect that the binding of the hidden charmed strange system \(D^*{{\bar{D}}}_{s}+D{\bar{D}}^*_s\) is loosely than the hidden charmed system \(D^{(*)}{{\bar{D}}}^{(*)}\) because more light exchanges are allowed for the latter state.

Table 4 The binding energies of some bound states with selected value of cutoff \(\varLambda \). The result for the \(D^{(*)}\bar{{D}_{s}}^{(*)}\) systems are listed in the 2–4th columns. The results for the \(D^{(*)}{\bar{D}}^{(*)}\) only with \(J/\psi \) exchange are listed in 5–7th columns. The cutoff \(\varLambda \) and binding energy \(E_B\) are in the units of GeV and MeV

In Table 3 we also present results for the systems \({D_{s}}{{\bar{D}}}^*_{s}\)/\({B_{s}}{{\bar{B}}}^*_{s}\) with different C parities. One can find that the bound states are produced from both \({D_{s}}{{\bar{D}}}^*_{s}\) and \({B_{s}}{{\bar{B}}}^*_{s}\) interactions with two C parties at small cutoffs. The two bound states of the \({D_{s}}{{\bar{D}}}^*_{s}\) system appears at cutoff of about 1.5 GeV, and the \({B_{s}}{{\bar{B}}}^*_{s}\) system is bound at a smaller cutoff of about 1 GeV. For doubly heavy systems \({D_{s}}{{D}}^*_{s}\)/\({B_{s}}{{B}}^*_{s}\), the spin parity is \(1^+\). A bound state can be found at a cutoff of about 3.4 GeV from the \(D_sD^*_s\) interaction and a bound state can be produced at a cutoff of about 2.8 GeV from the \(B_sB^*_s\) interaction.

In Table 5, the binding energies of the bound states with two vector mesons, \(D^*_{(s)}{\bar{D}}^*_s/B^*_{(s)}{\bar{B}}^*_s\), and \(D^*_{(s)}{D}^*_s/B^*_{(s)}{B}^*_s\), are presented. Three spin parities \(0^+\), \(1^+\), and \(2^+\) can be produced from S-wave interaction of two vector mesons. For the \(D^*_s{\bar{D}}^*_s\) system, a bound state with \(1^+\) appears at cutoff of 1.7 GeV. The \(D^*_s{\bar{D}}^*_s\) bound states with \(0^+\) and \(2^+\) are also produced with a relatively small cutoff of 1.4 and 2.0 GeV, respectively, which is consistent with the prediction of existence of a \(0^{++}\) molecular state from the \(D^*_s{\bar{D}}^*_s\) interaction [46,47,48,49,50]. The bound states are produced from the \(D^*_s{\bar{D}}^*_s\) interaction with three spins, but cutoffs larger than 4 GeV are required. In the bottom sector, the \(B^*{\bar{B}}^*_s\) interaction is still attractive for three spins, but not strong enough to produce bound state at a cutoff smaller than 5 GeV. Three bound states can be found in the \(B^*_s{\bar{B}}^*_s\) interaction for three spins at small cutoff as the \(D^*_s{\bar{D}}^*_s\) interactions. The cutoffs to produce \({D}_{s}^{(*)}{{\bar{D}}}^{(*)}_{s}\) and \({{B}_s}^{(*)}{{\bar{B}}}^{(*)}_{s}\) states are much smaller than these for \(D^{(*)}{{\bar{D}}}^{(*)}_{s}\) and \(B^{(*)}{{\bar{B}}}^{(*)}_{s}\) states because exchanges of light meson \(\eta \) and \(\phi \) take part in former interactions.

Table 5 The binding energies of the bound states from the interactions \(D^*_{(s)}{\bar{D}}^*_s/B^*_{(s)}{\bar{B}}^*_s\), and \(D^*_{(s)}{D}^*_s/B^*_{(s)}{B}^*_s\) with some selected values of cutoff \(\varLambda \). Other notations are the same as Table 2

In the doubly charmed sector, no bound state is found from the \(D^*D^*_s\) and \(D^*_sD^*_s\) interactions with spin parity \(0^+\) while bound states can be produced at cutoffs of about 3 GeV for spin parities \(1^+\) and \(2^+\). In the doubly bottom sector, the interactions become more attractive, one can find small cutoffs are required compared with these in the charmed sector, and a bound state can be found from the \(B^*B^*_s\) interaction with \(0^+\).

In our calculation we adopt the widely used assumption that no \(s{\bar{s}}\) component in \(\sigma \) meson. Hence, it can not be exchanged between the mesons considered in the current work. In Ref. [68] where the \(J/\psi \) exchange was not included, the \(\sigma \) exchange is proposed to play the most important factor to form a molecular state to interpret the \({Z_{cs}}(3985)^-\). As a complement, the influence of \(\sigma \) exchange in our calculation is listed in Table 6. Here, the Lagrangians in Eq. (8) for the vertices \({{{\mathcal {P}}}}{{{\mathcal {P}}}}\sigma \) and \({{{\mathcal {P}}}}^*{{{\mathcal {P}}}}^*\sigma \) are also applied to the heavy-strange meson. The binding energies of bound states with selected value of cutoff \(\varLambda \) are presented. One can find that there is no significant difference between cases with and without \(\sigma \) exchange. It suggests that if we include the \(J/\psi \) exchange, the contribution from the \(\sigma \) exchange may be smeared in our theoretical frame.

Table 6 The binding energies of some bound states with some selected values of cutoff \(\varLambda \). The results with \(\sigma \) exchange are listed 3–4th columns. The result without \(\sigma \) exchange are listed 5–6th columns. The cutoff \(\varLambda \) and binding energy \(E_B\) are in the units of GeV and MeV

3.2 Numerical results with coupled-channel calculation

Now, we consider the coupled-channel effect between the systems considered in the current work. The results for hidden charmed systems and doubly charmed systems are listed in Tables 7 and 8, respectively. Only the systems with the same quark constitutes and spin parity can be coupled. With couplings, the position z of states above the lowest threshold involved will acquire an imaginary part, which corresponds to the width as \(\varGamma =-2 \mathrm{Im} z\). To compare with the single-channel results, in the Tables we present the position as \(M_{th}-z\) instead of the position z of the pole, with the \(M_{th}\) being the nearest threshold.

Table 7 The \(M_{th}-z\) of the poles from the hidden-heavy coupled-channel interaction. The cutoff \(\varLambda \) and position of the pole z are in units of GeV and MeV, respectively

In the first part of Table 7, we present the results for \(D^{(*)}{\bar{D}}^{(*)}_s\) interactions. For the spin parity \(0^+\), only interactions \(D{\bar{D}}_s\) and \(D^*{\bar{D}}_s^*\) involve in the hidden charm strange system. In Tables 2 and 5, bound states are found at cutoff of 4.4 GeV and 4.2 GeV from the interactions \(D{\bar{D}}_s\) and \(D^*{\bar{D}}_s^*\), respectively. After coupling effect is included, two poles are found near the \(D{\bar{D}}_s\) and \(D^*{\bar{D}}_s^*\) thresholds at about 4.2 GeV. In the bottom sector, as in the single-channel calculation, no bound state is produced from the \(B{\bar{B}}_s-B^*{\bar{B}}_s^*\) interaction in the range of cutoff considered. For hidden charm hidden-strange system \(D_s{\bar{D}}_s-D^*_s{\bar{D}}^*_s\) with \(0^+\), two poles are produced near two thresholds at a cutoff of about 1.5 GeV, and the variation of the position of lower pole is much slower than the higher one, which is consistent with the results in the single-channel calculation. The result for the system \(B_s{\bar{B}}_s-B^*_s{\bar{B}}^*_s\) with \(0^+\) is analogous but with smaller cutoff.

For hidden charmed strange system with \(1^+\), three channels \(D^*{\bar{D}}_s\), \(D{\bar{D}}^*_s\), and \(D^*{\bar{D}}_s^*\) involve. In the single-channel calculation, we study the former two channels by constructing a U/V spin wave function under the SU(3)\(_F\) symmetry. There, the masses of \(D_s\) and D, as well as \(D^*_s\) and \(D^*\) mesons, will be chosen as their average values. However, these masses are different in experiment, which leads to larger violation of symmetry than the C parity in Eq. (3) where the difference between the masses of \(D^{(*)}_s\) mesons with different charges is very small in experiment. Here we take wave functions in Eq. (2) to perform the coupled-channel calculation to discuss the effect of violation of symmetry, and results for \(D^*{\bar{D}}_s-D{\bar{D}}^*_s-D^*{\bar{D}}^*_s\) system with \(1^+\) are given in Table 7. Generally speaking, the results with such treatment are consistent with single-channel calculation. Poles are produced near the \(D^*{\bar{D}}_s\) and \(D{\bar{D}}^*_s\) thresholds but with a little larger cutoff, which reflects effect of violation of the SU(3)\(_F\) symmetry. The pole near \(D^*{\bar{D}}_s^*\) threshold appears at a cutoff of about 4 GeV, which is analogous to the single-channel result.

For the hidden charm hidden strange systems, the wave functions can be constructed with fixed C parity, which can be well defined. Since the \(D_s^*{\bar{D}}_s^*\) has a \(C=-\) parity, we only consider \(D^*_s{\bar{D}}_s-D^*_s{\bar{D}}^*_s\) system with \(1^{+-}\). Two poles are reproduced near \(D^*_s{\bar{D}}_s\) and \(D^*_s{\bar{D}}^*_s\) thresholds at a cutoff of about 1.4 GeV. The above results suggest that inclusion of couplings between the channels considered in the current work affects the single-channel results very small. And violation of SU(3)\(_F\) symmetry provides visible effect, but the conclusion is unchanged.

In Table 8, the results for doubly charmed and doubly-bottomed systems with coupled-channel calculation are presented. No pole is found from all systems considered with \(0^+\), which is consistent with the single-channel results. For strange system \(DD_s^*-D^*D_s-D^*D_s^*\), poles are found near the \(D^*D_s\) and \(DD^*_s\) thresholds with spin parity \(1^+\) at cutoff of about 2.2 and 2.3 GeV, respectively. In doubly bottom sector, the poles are found at cutoffs of about 1.4 and 1.2 GeV, respectively. As in the hidden heavy sector (see Table 7), such results are different from the single-channel results with U/V spins quantitatively, but the conclusion is analogous qualitatively. The poles near the \(D^*D^*_s\) and \(B^*B^*_s\) thresholds appear at cutoffs of 2.5 and 1.2 GeV, which is similar to the values in the single-channel calculation. For the hidden-strange system, the poles are found near two thresholds of \(D^*_sD_s-D^*_sD^*_s\) system at cutoffs of 2.3 GeV and 3.3 GeV, and two thresholds of \(B^*_sB_s-B^*_sB^*_s\) system at cutoffs of 1.2 and 2.6 GeV, respectively.

Table 8 The \(M_{th}-z\) of the poles from the doubly heavy coupled-channel interaction. The cutoff \(\varLambda \) and \(M_{th}-z\) are in units of GeV and MeV, respectively

4 Summary and discussion

Inspired by the newly observed \({Z_{cs}}(3985)^-\), we study possible hidden and doubly heavy molecular states with hidden and open strangeness from interactions of \(D^{(*)}{{\bar{D}}}^{(*)}_{s}\)/\(B^{(*)}{{\bar{B}}}^{(*)}_{s}\), \({D}^{(*)}_{s}{{\bar{D}}}^{(*)}_{s}\)/\({{B}}^{(*)}_{s}{{\bar{B}}}^{(*)}_{s}\), \({D}^{(*)}{D_{s}}^{(*)}\)/\({B}^{(*)}{B_{s}}^{(*)}\) and \({D_{s}}^{(*)}{D_{s}}^{(*)}\)/ \({B_{s}}^{(*)}{B_{s}}^{(*)}\) in a qBSE approach. In the interactions of these systems, the light meson exchanges, as well as the \(J/\psi /\varUpsilon \) exchange, are included to construct the potential to find molecular state as a pole of the scattering amplitude.

The \(D{\bar{D}}^*_s\) interaction is considered with the U/V spins, a bound state is produced at a cutoff of about 3 GeV with spin parity \(1^+\) and \(G_{U/V}=+\), which can be related to the experimentally observed \(Z_{cs}(3985)\). An obvious difference between the hidden charmed interaction and hidden charm strange interaction is that only \(J/\psi \) exchange happens for the latter. If we take only the \(J/\psi \) exchange, an explicit calculation suggests that very similar binding energies can be obtained for two systems, which suggests that the \(Z_{cs}(3985)\) can be seen as a partner of \(Z_c(3900)\) state with \(I^G(J^{PC})=1^+(1^{+-})\) in the molecular state picture. Since light meson exchange is possible for hidden charm state, the \(Z_c(3900)\) state is more deeply bound than the \(Z_{cs}(3985)\) in molecular state picture. The interaction of \(D^*{\bar{D}}^*_s\) and its bottom partner are also attractive, but much weaker than \(D{\bar{D}}^*_s\) interaction, which makes it difficult to produce a bound state at a cutoff required to reproduce the the \(Z_{cs}(3985)\).

Besides above states, the calculation favors the existence of hidden-heavy states \(D_s{\bar{D}}_s/B_s{\bar{B}}_s\) with \(0^+\), \(D_s{\bar{D}}^*_s/B_s{\bar{B}}^*_s\) with \(1^{+\pm }\), \(D^*_s{\bar{D}}^*_s/B^*_s{\bar{B}}^*_s\) with \(0^+\), \(1^+\), and \(2^+\). The experimental search for the \(D^*_s{\bar{D}}^*_s/B^*_s{\bar{B}}^*_s\) molecular states is strongly suggested by our results. In the doubly heavy sector, the bound states can be found from the interactions \(DD^*_s/BB^*_s\) with \(1^+\) and \(G'=-\), \(D_s{\bar{D}}_s^*/B_s{\bar{B}}_s^*\) with \(1^+\), \(D^*D^*_s/B^*B^*_s\) with \(1^+\) and \(2^+\), and \(D^*_sD^*_s/B^*_sB^*_s\) with \(1^+\) and \(2^+\). Some other interactions are also found attractive, but may be not strong enough to produce a bound state.

In the calculation, we discuss the roles of the meson exchanges. The \(J/\psi /\varUpsilon \) exchange is found very important to form a molecular state. The contribution of such exchange overwhelms the contributions from light meson exchanges, as well as the \(\sigma \) exchange in our approach, which need future studies. In Ref. [27], an analysis also suggests that the \(J/\psi /\varUpsilon \) exchange is equivalent to the contact interaction in other approach [26]. The current results seem to show the importance of short-range interaction in formation of a molecular states, which is consistent with the conclusion in our previous work [27, 54] and other approach [26], where the short-range interaction is found essential to reproduce the \(Z_c(3900)\). We also perform a coupled-channel calculation to check the effect of coupling between different interactions. The results suggest that the conclusion is almost unchanged compared with the single-channel calculation. Of course, in the current work, only channels composed of open heavy flavor mesons are considered, further studies are required to give more precise prediction.