#  $Z_{cs}$, $Z_{c}$ and $Z_{b}$ states under the complex scaling method

Jian-Bo Cheng1<EMAIL_ADDRESS>Bo-Lin Huang1<EMAIL_ADDRESS>Zi-Yang
Lin1<EMAIL_ADDRESS>Shi-Lin Zhu1<EMAIL_ADDRESS>1School of Physics and
Center of High Energy Physics, Peking University, Beijing 100871, China

###### Abstract

We investigate the $Z_{b}$, $Z_{c}$ and $Z_{cs}$ states within the chiral
effective field theory framework and the $S$-wave single channel molecule
picture. With the complex scaling method, we accurately solve the Schrödinger
equation in momentum space. Our analysis reveals that the $Z_{b}(10610)$,
$Z_{b}(10650)$, $Z_{c}(3900)$ and $Z_{c}(4020)$ states are the resonances
composed of the $S-$wave $(B\bar{B}^{*}+B^{*}\bar{B})/\sqrt{2}$,
$B^{*}\bar{B}^{*}$, $(D\bar{D}^{*}+D^{*}\bar{D})/\sqrt{2}$ and
$D^{*}\bar{D}^{*}$, respectively. Furthermore, although the $Z_{cs}(3985)$ and
$Z_{cs}(4000)$ states exhibit a significant difference in width, these two
resonances may originate from the same channel, the $S-$wave
$(D_{s}\bar{D}^{*}+D_{s}^{*}\bar{D})/\sqrt{2}$. Additionally, we find two
resonances in the $S-$wave $D_{s}^{*}\bar{D}^{*}$ channel, corresponding to
the $Z_{cs}(4123)$ and $Z_{cs}(4220)$ states that await experimental
confirmation.

##  I Introduction

In the past decade, ongoing experimental efforts have led to the discovery of
a series of heavy quarkonium-like states known as the $XYZ$ states. The
charged $Z$ states like $Z_{c}(3900)$ and $Z_{c}(4020)$ provide strong
evidence of the exotic states, as they involve the light quarks to explain
their non-zero electric charge. Experimental advancements in the $Z_{b}$
sector can be traced back to 2011 when the Belle collaboration reported two
charged exotic candidates, $Z_{b}(10610)$ and $Z_{b}(10650)$ [1], which were
later confirmed in subsequent studies [2, 3]. Multiple hidden-charm tetraquark
candidates of the $Z_{c}$ states have been observed by the BESIII, Belle and
CLEO collaborations in electron-positron annihilation, including the charged
and neutral $Z_{c}(3900)$ and $Z_{c}(4020)$ states [4, 5, 6, 7, 8, 9, 10, 11,
12, 13, 14]. These states, with their masses near the thresholds of
$B^{(*)}\bar{B}^{*}$ and $D^{(*)}\bar{D}^{*}$, have been widely interpreted as
the molecule states in the papers [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25,
26, 27, 28, 29, 30]. Additionally, the existence of the strange partners with
the $Q\bar{Q}s\bar{q}^{\prime}$ ($q,q^{\prime}=u,d$) configurations is
predicted by the SU(3)-flavor symmetry, and indeed they have been discovered
in recent years.

In 2021, the BESIII collaboration observed an exotic hadron near the mass
thresholds of $D_{s}^{-}D^{*0}$ and $D_{s}^{*-}D^{0}$ in the processes
$e^{+}e^{-}\to K^{+}D_{s}^{-}D^{*0}$ and $K^{+}D_{s}^{*-}D^{0}$ [31]. The
corresponding mass and width fitted with a Breit-Wigner line shape are

$\displaystyle M[Z_{cs}(3985)]=3982.2_{-2.6}^{+1.8}\pm 2.1\ \text{MeV}\
\text{and}$ $\displaystyle\Gamma[Z_{cs}(3985)]=12.8_{-4.4}^{+5.3}\pm 3.0\
\text{MeV}.$ (1)

Last year, they observed a neutral $Z_{cs}(3985)^{0}$ in the processes
$e^{+}e^{-}\to K_{S}^{0}D_{s}^{+}D^{*-}$ and $K_{S}^{0}D_{s}^{*+}D^{-}$ [32].
The mass and width of the neutral $Z_{cs}(3985)^{0}$ have been determined to
be ($3992.2\pm 1.7\pm 1.6$) MeV and ($7.7^{+4.1}_{-3.8}\pm 4.3$) MeV,
respectively. Its mass, width and cross section are similar to those of the
charged $Z_{cs}(3985)^{+}$, which suggests that the neutral $Z_{cs}(3985)^{0}$
is the isospin partner of the $Z_{cs}(3985)^{+}$. Furthermore, in 2021, the
LHCb collaboration reported a series of distinct $Z_{cs}$ states. In the
hidden charm decay process $B^{+}\to J/\psi\phi K^{+}$, they observed two
$Z_{cs}$ states with $J^{P}=1^{+}$ [33]. One of these $Z_{cs}$ states is the
$Z_{cs}(4000)^{+}$, which is discovered with high significance. Its mass and
width are measured to be

$\displaystyle M[Z_{cs}(4000)]=4003\pm 6_{-14}^{+4}\ \text{MeV}\ \text{and}$
$\displaystyle\Gamma[Z_{cs}(4000)]=131\pm 15\pm 26\ \text{MeV},$ (2)

respectively. Additionally, the other $Z_{cs}$ state, $Z_{cs}(4220)^{+}$, has
a mass of $4216\pm 24_{-30}^{+43}$ MeV and a width of $233\pm 52_{-73}^{+97}$
MeV. The LHCb collaboration considers the $Z_{cs}(4000)^{+}$ and
$Z_{cs}(3985)^{+}$ to be distinct states due to their apparently different
widths, despite their close mass.

This discovery of the exotic $Z_{cs}$ hadrons inspired various theoretical
interpretations, including the compact tetraquark picture [34, 35, 36], the
molecule picture [37, 38, 39, 40, 41, 42, 43, 44], the mixing scheme [45, 46,
47, 48, 49] and the cusp effect [50]. When examining the BESIII and LHCb
observations of the $Z_{cs}$ states, some authors of the Refs. [41, 51, 52]
proposed that the $Z_{cs}(3985)$ and $Z_{cs}(4000)$ are the same entity,
whereas the Refs. [44, 48, 35, 39, 34, 47] considered them to be distinct
hadrons. Moreover, one can gain further insights from the comprehensive
reviews published in recent years [53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63,
64, 65, 66].

In Refs. [35, 39], the authors considered the $Z_{cs}(3985)$ and
$Z_{cs}(4000)$ as the SU(3)-flavor partners of $Z_{c}(3900)$, whose neutral
nonstrange members have opposite $C$ parity. The authors suggested the
$Z_{cs}(4000)/Z_{cs}(3985)$ is the pure molecular state composed of
$(|\bar{D}_{s}^{*}D\rangle+/-|\bar{D}_{s}D^{*}\rangle)/\sqrt{2}$. In addition,
they also predicted the existence of a molecule composed of
$\bar{D}_{s}^{*}D^{*}$ which may be confirmed by the BESIII in the subsequent
experiment [67]. However, the huge difference of their widths seems still hard
to interpret.

In this study, we employ the chiral effective field theory (ChEFT) to
investigate the properties of the $Z_{b}$, $Z_{c}$ and $Z_{cs}$ states in the
molecular picture. To explore the existence and relationships of the possible
resonances, we utilize the complex scaling method (CSM) [68, 69], which is a
powerful tool that provides a consistent treatment of the bound states and
resonances. We focus solely on the $S-$wave open-charm interaction, while
neglecting the possible contributions from the hidden charm. As illustrated in
our previous works [70, 71], we consider the cross diagram
$D\bar{D}^{*}\leftrightarrow D^{*}\bar{D}$ of the one-pion-exchange (OPE)
contribution. This contribution introduces a complex potential arising from
the three-body decay effect, which we take into account when investigating the
widths of the resonances.

This paper is organized as follows. In Sec. II, we introduce our framework
explicitly. In Sec. III, we present the effective Lagrangians and potentials.
In Sec. IV, we solve the complex scaled Schrödinger equation and give the
results of the $Z_{b}$, $Z_{c}$ and $Z_{cs}$. The last section V is a brief
summary.

##  II Framework

In this study, we consider the $Z_{b}$, $Z_{c}$ and $Z_{cs}$ states as the
molecular systems with the quantum numbers $I^{G}(J^{PC})=1^{+}(1^{+-})$,
$I^{G}(J^{PC})=1^{+}(1^{+-})$ and $I(J^{P})=1/2(1^{+})$, respectively. The
specific molecule systems under investigation are
$(B\bar{B}^{*}+B^{*}\bar{B})/\sqrt{2}$, $B^{*}\bar{B}^{*}$,
$(D\bar{D}^{*}+D^{*}\bar{D})/\sqrt{2}$, $D^{*}\bar{D}^{*}$,
$(D_{s}\bar{D}^{*}+D_{s}^{*}\bar{D})/\sqrt{2}$ and $D_{s}^{*}\bar{D}^{*}$.

In the earlier work [17], the $Z_{b}$ states were proposed as the bound states
of $\left[B\bar{B}^{*}+B^{*}\bar{B}\right]/\sqrt{2}$ and $B^{*}\bar{B}^{*}$.
The authors considered the $D$-wave channel and found that the $S-D$ wave
mixing effect could contribute significantly. Recent experiments [2, 3] have
uncovered additional evidence supporting the interpretation of the $Z_{b}$
states as the resonances. These findings show that the masses of the $Z_{b}$
states are higher than the threshold of the $B^{(*)}\bar{B}^{*}$ pairs, and
they can decay into the $B^{(*)}\bar{B}^{*}$ channel with the partial widths
in the range of tens of MeV. These findings strongly favor the resonance
interpretation over the bound state scenario.

The present CSM work confirms that the $D$-wave channel has a minimal impact
on the mass and width of the states. Therefore, we ignore the $D$-wave channel
in this work. On the other hand, we find that the coupled channel effect
between $(B\bar{B}^{*}+B^{*}\bar{B})/\sqrt{2}$ and $B^{*}\bar{B}^{*}$ is
negligible for the near threshold states. In addition, there are inelastic
channels in the final decay process, like $\Upsilon(nS)\pi$, that could be the
constituents of the $Z_{b}$ states as well. However, the couplings strength
between the $Z_{b}$ and the hidden-bottom channels is apparently smaller than
that between $Z_{b}$ and the open-bottom channels. Therefore, the influence of
the correction from the hidden-bottom channels should not be significant.
Furthermore, for the $Z_{c}$ and $Z_{cs}$ systems, we adopt the same
assumption that the inelastic hidden-heavy channels are not the primary
constituents. As a result, we focus on the simplest case, considering only the
$S-$wave open-heavy single channel.

The masses of the charmed meson and exchanged light mesons are collected in
Table 1. We take the isospin average masses to deal with the isospin
conserving process.

Mesons | Mass(MeV) | Mesons | Mass(MeV)
---|---|---|---
$D^{+}$ | 1869.66 | $B^{*}$ | 5324.70
$D^{0}$ | 1864.84 | $D_{s}^{+}$ | 1968.34
$D^{*+}$ | 2010.26 | $D_{s}^{*+}$ | 2112.2
$D^{*0}$ | 2006.85 | $\pi^{\pm}$ | 139.57
$B^{+}$ | 5279.34 | $\pi^{0}$ | 134.98
$B^{0}$ | 5279.65 |  |
Table 1: The masses of the charmed, bottomed and pion mesons, which are taken
from Ref. [72].

###  II.1 A brief discussion on the CSM

We first provide a brief overview of the CSM proposed by Aguilar, Balslev and
Combes in the 1970s [68, 69], commonly known as the ABC theorem. The CSM is a
powerful approach that allows for the treatment of resonances in a manner
similar to the bound states. The transformation of the radial coordinate $r$
and its conjugate momentum $k$ in the CSM are defined by:

$\displaystyle U(\theta)r=re^{i\theta},\qquad U(\theta)k=ke^{-i\theta}.$ (3)

After the complex scaling operation, the Schrödinger equation

$\displaystyle\frac{p^{2}}{2m}\phi_{l}(p)+\int\frac{p^{\prime
2}dp^{\prime}}{(2\pi)^{3}}V_{l,l^{\prime}}(p,p^{\prime})\phi_{l^{\prime}}(p^{\prime})=E\phi_{l}(p)$
(4)

in the momentum space becomes

$\displaystyle\frac{p^{2}e^{-2i\theta}}{2m}\tilde{\phi}_{l}(p)+\int\frac{p^{\prime
2}e^{-3i\theta}dp^{\prime}}{(2\pi)^{3}}V_{l,l^{\prime}}(pe^{-i\theta},p^{\prime}e^{-i\theta})\tilde{\phi}_{l^{\prime}}(p^{\prime})$
$\displaystyle=E\tilde{\phi}_{l}(p),$ (5)

with the normalization relation

$\displaystyle\frac{e^{-3i\theta}}{(2\pi)^{3}}\int_{0}^{\infty}\tilde{\phi}_{l}(p)^{2}p^{2}dp=1,$
(6)

where $l,l^{\prime}$ are the orbital angular momenta, and $p$ represents the
momentum in the center-of-mass frame. The potential $V_{l,l^{\prime}}$ after
partial wave decomposition can be expressed as

$\displaystyle V_{l,l^{\prime}}$ $\displaystyle=$ $\displaystyle\int
d\bm{\Omega}^{\prime}\int
d\bm{\Omega}\sum_{m_{l^{\prime}}=-l^{\prime}}^{l^{\prime}}\langle
l^{\prime},m_{l^{\prime}};s,m_{j}-m_{l^{\prime}}|j,m_{j}\rangle$ (7)
$\displaystyle\times$ $\displaystyle\sum_{m_{l}=-l}^{l}\langle
l,m_{l};s,m_{j}-m_{l}|j,m_{j}\rangle\mathcal{Y}_{l^{\prime},m_{l^{\prime}}}^{*}(\theta^{\prime},\phi^{\prime})$
$\displaystyle\times$ $\displaystyle\mathcal{Y}_{l,m_{l}}(\theta,\phi)\langle
s,m_{j}-m_{l^{\prime}}|\mathcal{V}|s,m_{j}-m_{l}\rangle,$

where $s$ and $j$ represent the total spin and total angular momentum of
systems, $m_{l}$ is the corresponding magnetic quantum number. The
$\mathcal{Y}_{l,m_{l}}(\theta,\phi)$ represents the spherical harmonics
associated with the angular coordinates $\theta$, $\phi$. The potential
operator $\mathcal{V}$ acts on the states $|s,m_{j}-m_{l^{\prime}}\rangle$ and
$|s,m_{j}-m_{l}\rangle$.

After performing the complex scaling operation, the resonance pole crosses the
branch cut into the first Riemann sheet when the rotation angle $\theta$
reaches a sufficiently large value, as depicted in Fig. 1. Consequently, the
wave functions of the resonances become square-integrable, similar to those of
the normalizable bound states. Further information on this technique can be
found in Refs. [73, 74].

Figure 1: The eigenvalue distribution of the complex scaled Schrödinger
equation for the two-body systems.

###  II.2 Analyticity of the OPE potentials for the $DD^{*}$ system

In our previous works [70, 71], we investigated the double-charm tetraquark
system using the CSM method. Notably, we found that the $D\bar{D}^{*}$ system
exhibits a unique characteristic where the zeroth component of the transferred
momentum of the exchanged pion exceeds the pion mass. This leads to an
imaginary part in the OPE potential. If a pole is obtained in this system, it
would correspond to an energy with an imaginary part, which can be interpreted
as its half-width. In the current study, we encounter this situation when
examining the OPE potential in the $(D\bar{D}^{*}+D^{*}\bar{D})/\sqrt{2}$
system with $1^{+}(1^{+-})$.

When considering the process $D\bar{D}^{*}\to D^{*}\bar{D}$, one can get the
OPE potential as follows

$\displaystyle
V_{\pi}\propto\frac{g^{2}}{2f_{\pi}^{2}}\frac{(\bm{\epsilon}^{*}\cdot\bm{q})(\bm{\epsilon}\cdot\bm{q})}{\bm{q}^{2}+m_{\pi}^{2}-q_{0}^{2}},$
(8)

where $q$ represents the transferred momentum of the pion, and $q_{0}$ is its
zeroth component. Since $q_{0}\approx m_{D^{*}}-m_{D}>m_{\pi}$, the poles of
the OPE potential are located on the real transferred momentum axis. However,
when performing the integral along the real $p^{\prime}$ axis in Eq. (4), we
encounter a numerical divergence. Fortunately, the CSM can resolve this
divergence issue without altering the analyticity of the OPE potential.
Through a complex scaling operation, the pole of the OPE potential is rotated
away from the real momentum axis in the momentum plane. As a result, the
integral along the real momentum axis bypasses the pole, effectively avoiding
divergence.

As shown in Fig. 2, we denote the total energy of the
$D\bar{D}^{*}/D^{*}\bar{D}$ system as $E$ and assume the $D$ meson to be on-
shell. In this case, the expression for $q_{0}$ is given by
$q_{0}=E-\sqrt{m_{D}^{2}+\bm{p}^{2}}-\sqrt{m_{\bar{D}}^{2}+\bm{p^{\prime}}^{2}}$.
With the heavy quark approximation, we neglect the kinetic energy contribution
to $q_{0}$ and introduce an energy shift $E\to E+m_{D}+m_{D^{*}}$. As a
result, we obtain $q_{0}=E+m_{D^{*}}-m_{D}$.

Figure 2: Three-body intermediate diagram in the process $D\bar{D}^{*}\to
D^{*}\bar{D}$. The total energy of the $D\bar{D}^{*}/D^{*}\bar{D}$ is $E$, and
the mesons which are cut by the red dashed line are on shell.

In other processes, the three-body effect vanishes, and we should consider the
different values of $q_{0}$ in the OPE potential. The specific values of
$q_{0}$ for each case are summarized in Table 2.

Process | $D\bar{D}^{*}\to D^{*}\bar{D}$ | $D^{*}\bar{D}^{*}\to D^{*}\bar{D}^{*}$ | $B\bar{B}^{*}\to B^{*}\bar{B}$
---|---|---|---
$q_{0}$ | $E+m_{D^{*}}-m_{D}$ | $0$ | $m_{B^{*}}-m_{B}$
Table 2: The $q_{0}$ is the zeroth component of the transferred momentum. $E$
is the total energy relative to the corresponding threshold. The other cases
not listed all give $q_{0}=0$.

##  III Lagrangians and Potentials

For the interaction of two heavy mesons, the chiral effective Lagrangians are
constructed based on the heavy quark symmetry and SU(3)-flavor symmetry. The
explicit expressions are given by

$\displaystyle\mathcal{L}$ $\displaystyle=$ $\displaystyle-i\langle
H^{(Q)}_{b}v\cdot(\delta_{ba}\partial+i\Gamma_{ba})\bar{H}^{(Q)}_{a}\rangle+g\langle
H^{(Q)}_{b}\mathbb{A}_{ba}^{\mu}\gamma_{\mu}\gamma_{5}\bar{H}_{a}^{(Q)}\rangle$
$\displaystyle-i\langle\bar{\tilde{H}}_{b}^{(\bar{Q})}v\cdot(\delta_{ba}\partial+i\Gamma_{ba})\tilde{H}_{a}^{(\bar{Q})}\rangle+g\langle\bar{\tilde{H}}_{b}^{(\bar{Q})}\mathbb{A}_{ba}^{\mu}\gamma_{\mu}\gamma_{5}\tilde{H}_{a}^{(\bar{Q})}\rangle$

where $H^{(Q)}$ is defined as

$\displaystyle H_{a}^{(Q)}$ $\displaystyle=$
$\displaystyle\frac{1+\not{v}}{2}\left[P_{a}^{*\mu}\gamma_{\mu}-P_{a}\gamma_{5}\right].$
(10)

And $\bar{H}_{a}^{(Q)}$, $\tilde{H}^{(\bar{Q})}$ and
$\bar{\tilde{H}}_{a}^{(\bar{Q})}$ are

$\displaystyle\bar{H}_{a}^{(Q)}$ $\displaystyle=$
$\displaystyle\gamma_{0}H_{a}^{(Q)\dagger}\gamma_{0}=\left[P_{a}^{*\dagger\mu}\gamma_{\mu}+P_{a}^{\dagger}\gamma_{5}\right]\frac{1+\not{v}}{2},$
$\displaystyle\tilde{H}^{(\bar{Q})}$ $\displaystyle=$
$\displaystyle\left[\tilde{P}_{a}^{*\mu}\gamma_{\mu}+\tilde{P}_{a}\gamma_{5}\right]\frac{1-\not{v}}{2}\quad\text{and}$
$\displaystyle\bar{\tilde{H}}_{a}^{(\bar{Q})}$ $\displaystyle=$
$\displaystyle\gamma_{0}\tilde{H}^{(\bar{Q})\dagger}\gamma_{0}=\frac{1-\not{v}}{2}\left[\tilde{P}_{a}^{*\dagger\mu}\gamma_{\mu}-\tilde{P}_{a}^{\dagger}\gamma_{5}\right],$
(11)

respectively, with $P_{a}^{(*)}=\left(D^{(*)0},D^{(*)+},D_{s}^{(*)+}\right)$
and
$\tilde{P}_{a}^{(*)}=\left(D^{(*)-},\bar{D}^{(*)0},\bar{D}_{s}^{(*)-}\right)$.

The light meson concerned parts are given that

$\displaystyle\mathbb{A}_{\mu}$ $\displaystyle=$
$\displaystyle\frac{i}{2}[\xi^{\dagger}(\partial_{\mu}\xi)+(\partial_{\mu}\xi)\xi^{\dagger}],\quad\Gamma_{\mu}=\frac{i}{2}[\xi^{\dagger}(\partial_{\mu}\xi)-(\partial_{\mu}\xi)\xi^{\dagger}],$
$\displaystyle\xi$ $\displaystyle=$
$\displaystyle\text{exp}[\frac{i\mathcal{M}}{f_{\pi}}]\quad\text{and}$ (12)
$\displaystyle\mathcal{M}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{6}}&\pi^{+}&K^{+}\\\
\pi^{-}&\frac{\pi^{0}}{\sqrt{2}}+\frac{\eta}{\sqrt{6}}&K^{0}\\\
K^{-}&\bar{K}^{0}&-\frac{2}{\sqrt{6}}\eta\end{pmatrix},$ (13)

where the pion decay constant $f_{\pi}$ is equal to 132 MeV. The coupling
constant associated with the $\pi$ exchange is $g=0.59$ [75].

The corresponding OPE potential in momentum space can be expressed as follows

$\displaystyle V^{{D\bar{D}^{*}}/{B\bar{B}^{*}}}$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{2f_{\pi}^{2}}\frac{(\bm{\epsilon}^{*}\cdot\bm{q})(\bm{\epsilon}\cdot\bm{q})}{\bm{q}^{2}+m_{\pi}^{2}-q_{0}^{2}},$
$\displaystyle V^{D^{*}\bar{D}^{*}/B^{*}\bar{B}^{*}}$ $\displaystyle=$
$\displaystyle-\frac{g^{2}}{2f_{\pi}^{2}}\frac{(\bm{\mathcal{T}}_{1}\cdot\bm{q})(\bm{\mathcal{T}}_{2}\cdot\bm{q})}{\bm{q}^{2}+m_{\pi}^{2}-q_{0}^{2}},$

where $\bm{\mathcal{T}}_{1}$ and $\bm{\mathcal{T}}_{2}$ represent the spin 1
operator with the forms
$\bm{\mathcal{T}}_{1}=-i\bm{\epsilon}_{3}^{\dagger}\times\bm{\epsilon}_{1}$
and
$\bm{\mathcal{T}}_{2}=-i\bm{\epsilon}_{4}^{\dagger}\times\bm{\epsilon}_{2}$.
Since we focus solely on the $S$-wave interactions, we can replace the above
spin-dependent operator with
$(\bm{\epsilon}^{*}\cdot\bm{q})(\bm{\epsilon}\cdot\bm{q})\to\frac{1}{3}\bm{q}^{2}$
and
$(\bm{T}_{1}\cdot\bm{q})(\bm{T}_{2}\cdot\bm{q})\to\frac{1}{3}\bm{q}^{2}\bm{T}_{1}\cdot\bm{T}_{2}$.

Regarding the contact term interaction, we adopt the form derived in Ref.
[28]. Upon performing the partial wave decomposition, one can obtain the
$S$-wave contact potential as

$\displaystyle\left[V_{ct}\right]_{l,l^{\prime}}$ $\displaystyle=$
$\displaystyle\tilde{C}_{s}+C_{s}(p^{2}+p^{\prime 2}),$

where $\tilde{C}_{s}$ and $C_{s}$ represent the partial wave low energy
constants (LECs). We restrict our analysis to the lowest-order interaction and
do not consider higher-order effects, such as the one-loop contribution.

To obtain the effective potentials, we introduce a Gaussian regulator to the
potentials as follows

$\displaystyle V_{l,l^{\prime}}=V_{l,l^{\prime}}\exp\left(-\frac{p^{\prime
2}}{\Lambda^{2}}-\frac{p^{2}}{\Lambda^{2}}\right),$ (15)

where $\Lambda$ is the cutoff parameter. The parameters $\Lambda$,
$\tilde{C}_{s}$ and $C_{s}$ can be adjusted while keeping the coupling
constants in the OPE potential fixed.

##  IV Numerical results

During the numerical calculation process, we discretize the Schrödinger Eq.
(4) in momentum space using the Gaussian quadrature approach. We approximate
the integral over the potential as a weighted sum over $N$ integration points
for $p=k_{j}$ ($j=1,N$):

$\displaystyle\int_{0}^{\infty}dp^{\prime}p^{\prime
2}V(p,p^{\prime})\phi(p^{\prime})\simeq\sum_{j=1}^{N}\omega_{j}p_{j}^{2}V(p,p_{j})\phi(p_{j}),$
$\displaystyle\frac{p^{2}}{2m}\phi(p)+\frac{1}{(2\pi)^{3}}\sum_{j=1}^{N}\omega_{j}p_{j}^{2}V(p,p_{j})\phi(p_{j})=E\phi(p),$
(16)

where $p_{j}$ and $\omega_{j}$ represent the Gaussian quadrature points and
weights, respectively. Furthermore, for clarity, we will omit the orbital
angular momentum subscript from this point onward. In Eq. (16), we have $N$
unknowns $\phi(k_{j})$ and an unknown $\phi(k)$. To avoid the need to
determine the entire function $\phi(k)$, we restrict the solution to the same
values of $k_{i}$ used to approximate the integral. This lead to $N$ coupled
linear equations:

$\displaystyle\frac{p_{i}^{2}}{2m}\phi(p_{i})+\frac{1}{(2\pi)^{3}}\sum_{j=1}^{N}\omega_{j}p_{j}^{2}V(p_{i},p_{j})\phi(p_{j})=E\phi(p_{i}).$
(17)

Therefore, the Schrödinger equation can be expressed in matrix form as

$\displaystyle[H][\phi]=E[\phi],$ (18)

with explicit matrices form

$\displaystyle\begin{pmatrix}\frac{p_{1}^{2}}{2m}+\frac{1}{(2\pi)^{3}}\omega_{1}p_{1}^{2}V(p_{1},p_{1})&\frac{1}{(2\pi)^{3}}\omega_{2}p_{2}^{2}V(p_{1},p_{2})&\cdots&\frac{1}{(2\pi)^{3}}\omega_{N}p_{N}^{2}V(p_{1},p_{N})\\\
\frac{1}{(2\pi)^{3}}\omega_{1}p_{1}^{2}V(p_{2},p_{1})&\frac{p_{2}^{2}}{2m}+\frac{1}{(2\pi)^{3}}\omega_{2}p_{2}^{2}V(p_{2},p_{2})&\cdots&\\\
\vdots&&&\\\
\frac{1}{(2\pi)^{3}}\omega_{1}p_{1}^{2}V(p_{N},p_{1})&\cdots&\cdots&\frac{p_{N}^{2}}{2m}+\frac{1}{(2\pi)^{3}}\omega_{N}p_{N}^{2}V(p_{N},p_{N})\end{pmatrix}\begin{pmatrix}\phi(p_{1})\\\
\phi(p_{2})\\\ \vdots\\\ \phi(p_{N})\\\
\end{pmatrix}=E\begin{pmatrix}\phi(p_{1})\\\ \phi(p_{2})\\\ \vdots\\\
\phi(p_{N})\\\ \end{pmatrix},$

where the wave function $\phi(k)$ on the grid can be represented as the
$N\times 1$ vector
$[\phi(p_{i})]=\begin{pmatrix}\phi(p_{1})&\phi(p_{2})&\cdots&\phi(p_{N})\\\
\end{pmatrix}^{T}$. Then, we can effectively solve Eq. (4). To find solutions
for the complex Schrödinger Eq. (5), we can simply make the substitutions
$p_{i}\to p_{i}e^{-i\theta}$, $\omega_{i}\to\omega_{i}e^{-i\theta}$ and
$\phi(p_{i})\to\tilde{\phi}(p_{i})$.

###  IV.1 The $Z_{b}$ and $Z_{c}$ system

Figure 3: The wave functions $\tilde{\phi}(p)$ of the $Z_{c}$ state with the
$I^{G}(J^{PC})=1^{+}(1^{+-})$. The rotation angle $\theta=35^{\circ}$ and the
parameters $\Lambda=0.300$ GeV, $\tilde{C}_{s}=2.86\times 10^{2}$ GeV-2 and
$C_{s}=-59.9\times 10^{2}$ GeV-4. The two diagrams correspond to system (a)
$\left[D\bar{D}^{*}+D^{*}\bar{D}\right]/\sqrt{2}$ and (b) $D^{*}\bar{D}^{*}$
respectively.

In this subsection, we investigate the exotic hadrons $Z_{b}$ and $Z_{c}$
using ChEFT. A similar study has been performed in Ref. [28], where the
$Z_{c}(3900)$ and $Z_{c}(4020)$ ($Z_{b}(10510)$ and $Z_{b}(10650)$) are
interpreted as $\left[D\bar{D}^{*}+D^{*}\bar{D}\right]/\sqrt{2}$ and
$D^{*}\bar{D}^{*}$ ($\left[B\bar{B}^{*}+B^{*}\bar{B}\right]/\sqrt{2}$ and
$B^{*}\bar{B}^{*}$) molecule with $J^{P}=1^{+}(1^{+-})$, respectively.
However, in our present work using CSM, we find that the contributions from
the $D$-wave constituents are negligible. Therefore, we neglect the $S-D$
mixing effect and solely focus on the $S$-wave channel in this section.

In our analysis, as shown in Table 3, we perform a fit of the LECs for the two
$Z_{b}$ and $Z_{c}$ states. Comparing our results with those in Ref. [28], we
find a similar cutoff value $\Lambda$ within a reasonable range. However, the
LECs $\tilde{C}_{s}$ and $C_{s}$ exhibit some variations, which could be
attributed to our omissions of the $D$-wave channel and the higher order
contribution. Additionally, we calculate the root-mean-square (RMS) radii, as
shown in Table 3, and find that the sizes of the $Z_{b}$ states are smaller
than those of the $Z_{c}$ states. Interestingly, the sizes of the two $Z_{b}$
($Z_{c}$) states are nearly identical. Moreover, the corresponding wave
functions, as depicted in Fig. 5, exhibit a striking resemblance. This
phenomenon is reasonable since our analysis in this work does not account for
the higher-order spin-dependent correction terms. The satisfaction of the
heavy quark spin symmetry justifies the similarities in the energy, decay
width, size and wave function observed in the $Z_{b}$ and $Z_{c}$ states.

As discussed in Ref. [70], the $DD^{*}/D\bar{D}^{*}$ system considered as the
$T_{cc}^{+}/X(3872)$ state can decay into the three-body open-charm channels
$DD\pi/D\bar{D}\pi$. In the case of the isovector $D\bar{D}^{*}$ system, it is
also necessary to consider the influence of the three-body decay. The
numerical results in the scheme we adopt, shown in Table 3 (row ”Adopt”), are
very close to the results under the instantaneous approximation $q_{0}=0$.
This implies that the mass, width and size have minimal changes. The reason
why the choice of $q_{0}$ matters for the $T_{cc}^{+}$ system but not for the
$Z_{c}$ system can be understood as follows. The mass of the $T_{cc}^{+}$
state is below the threshold of the $DD^{*}$ system, making the two-body decay
process kinetically forbidden. Therefore, the three-body decay becomes the
dominant decay modes, and the value of $q_{0}$, which partly reflects the
three-body decay width, becomes important. On the other hand, the
$Z_{c}(3900)^{+}$ state is clearly above the threshold of the $D\bar{D}^{*}$
system, allowing for the two-body decay process. Since the contribution from
the three-body decay is significantly smaller in this case, the choice of
$q_{0}$ does not significantly alter the results.

System | Threshold |  $\left[m,\Gamma\right]_{\text{pole}}$(MeV) |  $\left[m,\Gamma\right]_{\text{exp}}$(MeV) | RMS(fm)
---|---|---|---|---
$\left[B\bar{B}^{*}+B^{*}\bar{B}\right]/\sqrt{2}$ | 10604.2 | $\left[10606.9_{-1.5}^{+1.8},15.0_{-3.2}^{+3.4}\right]$ | $\left[10607.2_{-2.0}^{+2.0},18.4_{-2.4}^{+2.4}\right]$ | $0.70_{-0.01}^{+0.07}-0.15_{-0.10}^{+0.09}i$
$B^{*}\bar{B}^{*}$ | 10649.4 | $\left[10652.2_{-1.6}^{+1.8},14.8_{-3.2}^{+3.4}\right]$ | $\left[10652.2_{-1.5}^{+1.5},11.5_{-2.2}^{+2.2}\right]$ | $0.70_{-0.02}^{+0.07}-0.15_{-0.11}^{+0.09}i$
$\left[D\bar{D}^{*}+D^{*}\bar{D}\right]/\sqrt{2}$(Adopt) | 3875.8 | $\left[3884.3_{-0.6}^{+0.6},26.0_{-1.4}^{+1.4}\right]$ | $\left[3881.7_{-2.3}^{+2.3},26.6_{-3.4}^{+3.0}\right]$ | $1.21_{-0.05}^{+0.06}+0.12_{-0.03}^{+0.03}i$
$\left[D\bar{D}^{*}+D^{*}\bar{D}\right]/\sqrt{2}$(Inst) | 3875.8 | $\left[3884.8_{-0.6}^{+0.6},25.8_{-1.4}^{+1.4}\right]$ | $\left[3881.7_{-2.3}^{+2.3},26.6_{-3.4}^{+3.0}\right]$ | $1.20_{-0.05}^{+0.06}+0.13_{-0.03}^{+0.03}i$
$D^{*}\bar{D}^{*}$ | 4017.1 | $\left[4025.8_{-0.6}^{+0.6},24.0_{-1.4}^{+1.3}\right]$ | $\left[4025.5_{-5.6}^{+3.7},26.0_{-6.0}^{+6.0}\right]$ | $1.20_{-0.05}^{+0.06}+0.13_{-0.03}^{+0.03}i$
Table 3: The extracted poles for all states are listed with the quantum numbers $I^{G}(J^{PC})=1^{+}(1^{+-})$. The fitted parameter for the $B^{*}\bar{B}^{(*)}$ system are $\Lambda=0.510^{+0.027}_{-0.041}$ GeV, $\tilde{C}_{s}=0.48^{+0.15}_{-0.13}\times 10^{2}$ GeV-2 and $C_{s}=-5.4^{+0.63}_{-0.65}\times 10^{2}$ GeV-4. The fitted parameter for the $D^{*}\bar{D}^{(*)}$ system are $\Lambda=0.300^{+0.012}_{-0.013}$ GeV, $\tilde{C}_{s}=2.86^{+0.21}_{-0.22}\times 10^{2}$ GeV-2 and $C_{s}=-59.9^{+2.8}_{-3.1}\times 10^{2}$ GeV-4. The RMS is the root-mean-square radius in the CSM, which has been discussed in the Ref. [76]. Its real part is interpreted as an expectation value, and the imaginary part corresponds to a measure of the uncertainty in observation. The data of row Adopt are the results we actually adopt, the $q_{0}$ herein is from Table 2. The data of row Inst are from the instantaneous approximation $q_{0}=0$. System | Threshold |  $\left[m,\Gamma\right]_{\text{pole}}$(MeV) |  $\left[m,\Gamma\right]_{\text{exp}}$(MeV) | RMS(fm)
---|---|---|---|---
$\left[D_{s}\bar{D}^{*}+D_{s}^{*}\bar{D}\right]/\sqrt{2}$(1*) | 3976.1 | $\left[3982.4^{+2.2}_{-2.1},14.1^{+3.7}_{-3.6}\right]$ | $\left[3982.5_{-3.3}^{+2.8},12.8_{-5.3}^{+6.1}\right]$ | $1.89^{+0.13}_{-0.11}+0.43^{+0.09}_{-0.14}i$
$\left[D_{s}\bar{D}^{*}+D_{s}^{*}\bar{D}\right]/\sqrt{2}$(2*) |  | $\left[4010.7^{+6.3}_{-6.2},119.6^{+14.5}_{-14.7}\right]$ | $\left[4003_{-15.2}^{+7.2},131_{-30.0}^{+30.0}\right]$ | $1.78^{+0.18}_{-0.14}+1.31^{+0.08}_{-0.07}i$
$D_{s}^{*}\bar{D}^{*}$(1) | 4119.1 | $\left[4125.2_{-2.1}^{+2.2},13.2_{-3.4}^{+3.5}\right]$ | $\left[4123.5_{-1.3}^{+1.3},-\right]$ | $1.89_{-0.11}^{+0.13}+0.43_{-0.14}^{+0.09}i$
$D_{s}^{*}\bar{D}^{*}$(2) |  | $\left[4152.7_{-6.0}^{+6.1},115.0_{-14.2}^{+14.0}\right]$ | $\left[4216_{-38}^{+49},233_{-90}^{+110}\right]$ | $1.78_{-0.14}^{+0.18}+1.31_{-0.07}^{+0.08}i$
Table 4: The poles are all listed with the number $I(J^{P})=1/2(1^{+})$. The
fitted parameter for the $D_{s}^{*}\bar{D}^{(*)}$ system are
$\Lambda=0.192^{+0.012}_{-0.013}$ GeV, $\tilde{C}_{s}=6.8^{+2.8}_{-2.7}\times
10^{2}$ GeV-2 and $C_{s}=-186.9^{+50.4}_{-64.4}\times 10^{2}$ GeV-4. The RMS
is the root-mean-square radius in the CSM, which has been discussed in the
Ref. [76]. Its real part is interpreted as an expectation value, and the
imaginary part corresponds to a measure of the uncertainty in observation. The
states labeled as “1*” and “2*” correspond to the input states. The symbol “-”
indicates that the width of the $Z_{cs}(4123)$ state has not been confirmed by
experiment yet.

###  IV.2 The $Z_{cs}$ system

Figure 4: The eigenvalue distribution of the $Z_{cs}$ with the
$I(J^{P})=1/2(1^{+})$. The parameters $\Lambda=0.192$ GeV,
$\tilde{C}_{s}=6.8\times 10^{2}$ GeV-2 and $C_{s}=-186.9\times 10^{2}$
GeV-4.The orange (green) points (square points) and lines correspond to the
complex rotation angle $\theta=35^{\circ}$ ($40^{\circ}$). The two diagrams
correspond to system: (a)
$\left[D_{s}\bar{D}^{*}+D_{s}^{*}\bar{D}\right]/\sqrt{2}$ (b)
$D_{s}^{*}\bar{D}^{*}$.

In Refs. [35, 39], the $Z_{cs}(3985)$ and $Z_{cs}(4000)$ states are discussed
as the SU(3)-flavor partners of $Z_{c}(3900)$, with their neutral nonstrange
members having opposite C parity. The authors suggest that the
$Z_{cs}(4000)/Z_{cs}(3985)$ state can be described as a pure molecular state
composed of $(|D_{s}\bar{D}^{*}\rangle+/-|D_{s}^{*}\bar{D}\rangle)/\sqrt{2}$.
Furthermore, they also predicted the existence of a $D_{s}^{*}\bar{D}^{*}$
molecular state, which is potentially supported by the recent work of the
BESIII collaboration [67]. These studies provide interesting insights into the
nature and composition of the $Z_{cs}$ states under the molecule picture.

However, an issue that remains unresolved is the significant difference in the
widths between the $Z_{cs}(3985)$ and the $Z_{cs}(4000)$. To address this
difference, we propose an alternative explanation where these two states are
considered as two resonances associated with the same system, namely
$(D_{s}\bar{D}^{*}+D_{s}^{*}\bar{D})/\sqrt{2}$. According to our proposal, the
$Z_{cs}(3985)$ corresponds to the resonance with a narrower width, while the
$Z_{cs}(4000)$ corresponds to the resonance with a broader width, as
illustrated in Fig. 4. This interpretation differs from the prevailing
viewpoints in literature.

Figure 5: The wave functions $\tilde{\phi}(p)$ of the $Z_{c}s$ state with the
$I(J^{P})=1/2(1^{+})$. The rotation angle $\theta=40^{\circ}$ and the
parameters $\Lambda=0.192$ GeV, $\tilde{C}_{s}=6.8\times 10^{2}$ GeV-2 and
$C_{s}=-186.9\times 10^{2}$ GeV-4. The four diagrams correspond to: (a) pole 1
of $\left[D_{s}\bar{D}^{*}+D_{s}^{*}\bar{D}\right]/\sqrt{2}$ system (b) pole 2
of $\left[D_{s}\bar{D}^{*}+D_{s}^{*}\bar{D}\right]/\sqrt{2}$ system (c) pole 1
of $D_{s}^{*}\bar{D}^{*}$ system (d) pole 2 of $D_{s}^{*}\bar{D}^{*}$ system.

In the previous subsection IV.1, we found that the
$\left[D\bar{D}^{*}+D^{*}\bar{D}\right]/\sqrt{2}$ and $D^{*}\bar{D}^{*}$
systems, associated with the two $Z_{c}$ states, exhibit similar outcomes due
to the heavy quark spin symmetry. Thus, it is feasible to employ the same
parameters for them. Following this scheme, we use the same parameters for
both the $(D_{s}\bar{D}^{*}+D_{s}^{*}\bar{D})/\sqrt{2}$ and
$D_{s}^{*}\bar{D}^{*}$ systems. By adopting the available experimental data of
$Z_{cs}(3985)$ and $Z_{cs}(4000)$, we determine the central values and errors
of $\Lambda$, $\tilde{C}{s}$ and $C_{s}$, and perform calculations for the
$D_{s}^{*}\bar{D}^{*}$ system. The corresponding parameter values, masses,
widths and sizes can be found in Table 4.

In the framework of ChEFT, it is generally expected that the cutoff region
should exceed the pion mass $m_{\pi}$ while not significantly exceeding $0.5$
GeV, as the higher-mass mesons ($\sigma$, $\rho$, $\omega$, etc.) are
integrated out. Consequently, the $\Lambda$ adopted in this study, $0.3\sim
0.5$ GeV, is reasonable for the $Z_{c}$ and $Z_{b}$ cases. However, in the
case of $Z_{cs}$, the OZI suppression prohibits the contributions from either
OPE or one-kaon-exchange. As a result, the contact term becomes the only
interaction that needs to be considered. This can be viewed as effectively
integrating out the pion and kaon fields. Therefore, we adopt a smaller value
of $\Lambda\approx 0.2$ GeV for the $Z_{cs}$ cases.

According to the results in Table 4, the newly reported $Z_{cs}(4123)$ by
BESIII collaboration [67] could correspond to the narrower
$D_{s}^{*}\bar{D}^{*}$ state, although the experimental width is yet to be
confirmed. Its estimated mass is around 4125.2 MeV and width is approximately
13.2 MeV. Furthermore, the $Z_{cs}(4220)$ is anticipated to correspond to a
broader resonance with its central values of the mass and width at 4152.7 MeV
and 115.0 MeV. Indeed, the mass and width of $Z_{cs}(4220)$ both fall within
the two-standard-deviation region of the experimental data.

Furthermore, as shown in Fig. 5 and Table 4, the narrow (or broader)
resonances exhibit remarkably similar wave functions and sizes.

##  V Summary

In this study, we employ the ChEFT to investigate the hidden-heavy tetraquark
states with $I^{G}(J^{PC})=1^{+}(1^{+-})$ and the hidden-charm states with a
strange quark with $I(J^{P})=1/2(1^{+})$ in the molecule picture. The
couplings between the $S$-wave open-heavy channel and other channels, such as
the $D$-wave channel, the $S$-wave channel with different constituents, and
the hidden-heavy channels, are expected to be small. Therefore, we focus on
the $S-$wave open-heavy single channels:
$\left[D\bar{D}^{*}+D^{*}\bar{D}\right]/\sqrt{2}$, $D^{*}\bar{D}^{*}$,
$\left[B\bar{B}^{*}+B^{*}\bar{B}\right]/\sqrt{2}$, $B^{*}\bar{B}^{*}$,
$(D_{s}\bar{D}^{*}+D_{s}^{*}\bar{D})/\sqrt{2}$ and $D^{*}_{s}\bar{D}^{*}$.

We employ the effective Lagrangians based on heavy quark symmetry and chiral
symmetry, considering both contact and OPE diagrams. To investigate the
possible resonances, we adopt the CSM to consistently analyze the bound states
and resonances. In contrast to our previous work [70, 77], we perform the
momentum space Schrödinger equation and discretize it using the Gaussian
quadrature approach.

In our investigation of the $Z_{b}$ system, we fit experimental data to
extract resonance parameters within the molecule picture. With $\Lambda=0.510$
GeV, $\tilde{C}_{s}=0.48\times 10^{2}$ GeV-2 and $C_{s}=-5.4\times 10^{2}$
GeV-4, we obtain the mass and width values of 10606.9 MeV and 15.0 MeV for the
$\left[B\bar{B}^{*}+B^{*}\bar{B}\right]/\sqrt{2}$ resonance, while 10652.2 MeV
and 14.8 MeV for the $B^{*}\bar{B}^{*}$ resonance. The RMS radii for these two
resonances are both approximately $0.70-0.15i$ fm. Similarly, we perform
calculations for the $Z_{c}$ system in the $S-$wave $1^{+}(1^{+-})$ channels:
$\left[D\bar{D}^{*}+D^{*}\bar{D}\right]/\sqrt{2}$ and $D^{*}\bar{D}^{*}$.
Taking $\Lambda=0.300$ GeV, $\tilde{C}_{s}=2.86\times 10^{2}$ GeV-2 and
$C_{s}=-59.9\times 10^{2}$ GeV-4, we obtain the mass and width values of
3884.3 MeV and 26.0 MeV for the former resonance, while 4025.8 MeV and 24.0
MeV for the latter resonance. The RMS radii for both resonances are around
$1.20+0.13i$ fm. For the isovector
$\left[D\bar{D}^{*}+D^{*}\bar{D}\right]/\sqrt{2}$ system, we also consider the
influence of the three-body decay. However, the numerical results under the
instantaneous approximation with $q_{0}=0$ (as shown in the “Inst” row of
Table 3) are very close to the results of the “Adopt” row, indicating the
minimal changes in the mass, width and size. Thus, we conclude that the 2-body
decay process dominates the width of this resonance.

We consider the hidden-charm tetraquark states with a strange quark and
propose that the $Z_{cs}(3985)$ and $Z_{cs}(4000)$ resonances correspond to
the same channel $(D_{s}\bar{D}^{*}+D_{s}^{*}\bar{D})/\sqrt{2}$. Taking the
data of $Z_{cs}(3985)$ and $Z_{cs}(4000)$ as input, we extract the central
values and errors of the parameters $\Lambda$, $\tilde{C}_{s}$ and $C_{s}$.
With $\Lambda=0.192$ GeV, $\tilde{C}_{s}=6.8\times 10^{2}$ GeV-2 and
$C_{s}=-186.9\times 10^{2}$ GeV-4, we obtain the mass and width values of
3982.4 MeV and 14.1 MeV for the $Z_{cs}(3985)$, while 4010.7 MeV and 119.6 MeV
for the $Z_{cs}(4000)$. The corresponding RMS radii are $1.89+0.43i$ fm and
$1.78+1.31i$ fm, respectively. For the $D^{*}_{s}\bar{D}^{*}$ system, we adopt
the same parameters based on the heavy quark spin symmetry and also find two
resonances. The narrower resonance has a mass of 4125.2 MeV and a width of
13.2 MeV, which nicely matches the observed $Z_{cs}(4123)$ reported by the
BESIII collaboration [67]. Hence, we interpret it as the $Z_{cs}(4123)$,
although the experimental width is yet to be confirmed. On the other hand, the
broader resonance has a mass of 4152.7 MeV and a width of 115.0 MeV. We
interpret it as the $Z_{cs}(4220)$ observed by the LHCb collaboration [33], as
its mass and width fall within the two-standard-deviation region of the
experimental data.

In summary, we apply the ChEFT to investigate the $Z_{b}$, $Z_{c}$ and
$Z_{cs}$ states. Our analysis suggests that the $Z_{b}(10610)$,
$Z_{b}(10650)$, $Z_{c}(3900)$ and $Z_{c}(4020)$ can be interpreted as the
molecular states formed by the $S-$wave $B\bar{B}^{*}$, $B^{*}\bar{B}^{*}$,
$D\bar{D}^{*}$ and $D^{*}\bar{D}^{*}$ constituents, respectively. Although the
$Z_{cs}(3985)$ and $Z_{cs}(4000)$ states exhibit a significant width
difference, these two resonances may originate from the same $S-$wave channel
$(D_{s}\bar{D}^{*}+D_{s}^{*}\bar{D})/\sqrt{2}$. We also find two resonances in
the $D_{s}^{*}\bar{D}^{*}$ channel, which can be identified as the
$Z_{cs}(4123)$ and $Z_{cs}(4220)$. Our results provide a prediction for the
width of the $Z_{cs}(4123)$ that awaits experimental confirmation.
Additionally, we offer a precise mass and width range for the $Z_{cs}(4220)$,
which can guide future experimental searches for the hidden-charm tetraquarks
with a strange quark.

###### Acknowledgements.

This research is supported by the National Science Foundation of China under
Grants No. 11975033, No. 12070131001 and No. 12147168. The authors thank K.
Chen, Y. Ma and B. Wang for helpful discussions.

## References

  * A. Bondar _et al_. (2012) [Belle Collaboration] A. Bondar _et al_. (Belle Collaboration), Phys. Rev. Lett. 108, 122001 (2012).
  * A. Garmash _et al_. (2015) [Belle Collaboration] A. Garmash _et al_. (Belle Collaboration), Phys. Rev. D 91, 072003 (2015).
  * A. Garmash _et al_. (2016) [Belle Collaboration] A. Garmash _et al_. (Belle Collaboration), Phys. Rev. Lett. 116, 212001 (2016).
  * M. Ablikim _et al_. (2013a) [BESIII Collaboration] M. Ablikim _et al_. (BESIII Collaboration), Phys. Rev. Lett. 110, 252001 (2013a).
  * M. Ablikim _et al_. (2013b) [BESIII Collaboration] M. Ablikim _et al_. (BESIII Collaboration), Phys. Rev. Lett. 111, 242001 (2013b).
  * M. Ablikim _et al_. (2014a) [BESIII Collaboration] M. Ablikim _et al_. (BESIII Collaboration), Phys. Rev. Lett. 112, 022001 (2014a).
  * M. Ablikim _et al_. (2014b) [BESIII Collaboration] M. Ablikim _et al_. (BESIII Collaboration), Phys. Rev. Lett. 112, 132001 (2014b).
  * M. Ablikim _et al_. (2014c) [BESIII Collaboration] M. Ablikim _et al_. (BESIII Collaboration), Phys. Rev. Lett. 113, 212002 (2014c).
  * M. Ablikim _et al_. (2015a) [BESIII Collaboration] M. Ablikim _et al_. (BESIII Collaboration), Phys. Rev. Lett. 115, 112003 (2015a).
  * M. Ablikim _et al_. (2015b) [BESIII Collaboration] M. Ablikim _et al_. (BESIII Collaboration), Phys. Rev. Lett. 115, 182002 (2015b).
  * M. Ablikim _et al_. (2015c) [BESIII Collaboration] M. Ablikim _et al_. (BESIII Collaboration), Phys. Rev. D 92, 092006 (2015c).
  * M. Ablikim _et al_. (2015d) [BESIII Collaboration] M. Ablikim _et al_. (BESIII Collaboration), Phys. Rev. Lett. 115, 222002 (2015d).
  * Z. Q. Liu _et al_. (2013) [Belle Collaboration] Z. Q. Liu _et al_. (Belle Collaboration), Phys. Rev. Lett. 110, 252002 (2013).
  * Xiao  _et al._ [2013] T. Xiao, S. Dobbs, A. Tomaradze, and K. K. Seth, Phys. Lett. B 727, 366 (2013).
  * Molina and Oset [2009] R. Molina and E. Oset, Phys. Rev. D 80, 114013 (2009).
  * Zhang  _et al._ [2011] J.-R. Zhang, M. Zhong, and M.-Q. Huang, Phys. Lett. B 704, 312 (2011).
  * Sun  _et al._ [2011] Z.-F. Sun, J. He, X. Liu, Z.-G. Luo, and S.-L. Zhu, Phys. Rev. D 84, 054002 (2011).
  * Ozpineci  _et al._ [2013] A. Ozpineci, C. W. Xiao, and E. Oset, Phys. Rev. D 88, 034018 (2013).
  * He  _et al._ [2013] J. He, X. Liu, Z.-F. Sun, and S.-L. Zhu, Eur. Phys. J. C 73, 2635 (2013).
  * Guo  _et al._ [2013] F.-K. Guo, C. Hidalgo-Duque, J. Nieves, and M. P. Valderrama, Phys. Rev. D 88, 054007 (2013).
  * Dong  _et al._ [2013] Y. Dong, A. Faessler, T. Gutsche, and V. E. Lyubovitskij, Phys. Rev. D 88, 014030 (2013).
  * Wang and Huang [2014] Z.-G. Wang and T. Huang, Eur. Phys. J. C 74, 2891 (2014).
  * He [2014] J. He, Phys. Rev. D 90, 076008 (2014).
  * Aceti  _et al._ [2014] F. Aceti, M. Bayar, E. Oset, A. M. Torres, K. P. Khemchandani, J. M. Dias, F. S. Navarra, and M. Nielsen, Phys. Rev. D 90, 016003 (2014).
  * Wang [2014] Z.-G. Wang, Eur. Phys. J. C 74, 2963 (2014).
  * Karliner and Rosner [2015] M. Karliner and J. L. Rosner, Phys. Rev. Lett. 115, 122001 (2015).
  * Baru  _et al._ [2019] V. Baru, A. A. Filin, C. Hanhart, A. V. Nefediev, Q. Wang, and J.-L. Wynen, EPJ Web Conf. 202, 06012 (2019).
  * Wang  _et al._ [2020] B. Wang, L. Meng, and S.-L. Zhu, Phys. Rev. D 102, 114019 (2020).
  * Ding  _et al._ [2020] Z.-M. Ding, H.-Y. Jiang, and J. He, Eur. Phys. J. C 80, 1179 (2020).
  * Dai  _et al._ [2022] L. R. Dai, E. Oset, A. Feijoo, R. Molina, L. Roca, A. M. Torres, and K. P. Khemchandani, Phys. Rev. D 105, 074017 (2022).
  * M. Ablikim _et al_. (2021) [BESIII Collaboration] M. Ablikim _et al_. (BESIII Collaboration), Phys. Rev. Lett. 126, 102001 (2021).
  * M. Ablikim _et al_. (2022) [BESIII Collaboration] M. Ablikim _et al_. (BESIII Collaboration), Phys. Rev. Lett. 129, 112003 (2022).
  * R. Aaij _et al_. (2021) [LHCb Collaboration] R. Aaij _et al_. (LHCb Collaboration), Phys. Rev. Lett. 127, 082001 (2021).
  * Wang  _et al._ [2023] Y.-H. Wang, J. Wei, C.-S. An, and C.-R. Deng, Chin. Phys. Lett. 40, 021201 (2023).
  * Maiani  _et al._ [2021] L. Maiani, A. D. Polosa, and V. Riquer, Science Bulletin 66, 1616 (2021).
  * Shi  _et al._ [2021] P.-P. Shi, F. Huang, and W.-L. Wang, Phys. Rev. D 103, 094038 (2021).
  * Chen and Huang [2021] R. Chen and Q. Huang, Phys. Rev. D 103, 034008 (2021).
  * Yan  _et al._ [2021] M.-J. Yan, F.-Z. Peng, M. Sánchez Sánchez, and M. Pavon Valderrama, Phys. Rev. D 104, 114025 (2021).
  * Meng  _et al._ [2021] L. Meng, B. Wang, G.-J. Wang, and S.-L. Zhu, Sci. Bull. 66, 2065 (2021).
  * Zhai  _et al._ [2022] Q.-Y. Zhai, M.-Z. Liu, J.-X. Lu, and L.-S. Geng, Phys. Rev. D 106, 034026 (2022).
  * Wu  _et al._ [2022] Q. Wu, D.-Y. Chen, W.-H. Qin, and G. Li, Eur. Phys. J. C 82, 520 (2022).
  * Meng  _et al._ [2022] L. Meng, B. Wang, G.-J. Wang, and S.-L. Zhu, Nuclear and Particle Physics Proceedings 318–323, 85 (2022).
  * Du  _et al._ [2022] M.-L. Du, M. Albaladejo, F.-K. Guo, and J. Nieves, Phys. Rev. D 105, 074018 (2022).
  * Chen [2022] H.-X. Chen, Phys. Rev. D 105, 094003 (2022).
  * Karliner and Rosner [2021] M. Karliner and J. L. Rosner, Phys. Rev. D 104, 034033 (2021).
  * Jin  _et al._ [2021] X. Jin, Y. Wu, X. Liu, H. Huang, J. Ping, and B. Zhong, Eur. Phys. J. C 81, 1108 (2021).
  * Yang  _et al._ [2021] G. Yang, J. Ping, and J. Segovia, Phys. Rev. D 104, 094035 (2021).
  * Han and Xiao [2022] S. Han and L.-Y. Xiao, Phys. Rev. D 105, 054008 (2022).
  * Cao  _et al._ [2023] Z.-H. Cao, W. He, and Z.-F. Sun, Phys. Rev. D 107, 014017 (2023).
  * Luo and Nakamura [2023] X. Luo and S. X. Nakamura, Phys. Rev. D 107, L011504 (2023).
  * Ortega  _et al._ [2021] P. G. Ortega, D. R. Entem, and F. Fernández, Phys. Lett. B 818, 136382 (2021).
  * Giron  _et al._ [2021] J. F. Giron, R. F. Lebed, and S. R. Martinez, Phys. Rev. D 104, 054001 (2021).
  * Hosaka  _et al._ [2016] A. Hosaka, T. Iijima, K. Miyabayashi, Y. Sakai, and S. Yasui, PTEP 2016, 062C01 (2016).
  * Chen  _et al._ [2016] H.-X. Chen, W. Chen, X. Liu, and S.-L. Zhu, Phys. Rep. 639, 1 (2016).
  * Chen  _et al._ [2017] H.-X. Chen, W. Chen, X. Liu, Y.-R. Liu, and S.-L. Zhu, Rept. Prog. Phys. 80, 076201 (2017).
  * Ali  _et al._ [2017] A. Ali, J. S. Lange, and S. Stone, Prog. Part. Nucl. Phys. 97, 123 (2017).
  * Esposito  _et al._ [2017] A. Esposito, A. Pilloni, and A. Polosa, Phys. Rep. 668, 1 (2017).
  * Lebed  _et al._ [2017] R. F. Lebed, R. E. Mitchell, and E. S. Swanson, Prog. Part. Nucl. Phys. 93, 143 (2017).
  * Guo  _et al._ [2018] F.-K. Guo, C. Hanhart, U.-G. Meißner, Q. Wang, Q. Zhao, and B.-S. Zou, Rev. Mod. Phys. 90, 015004 (2018).
  * Liu  _et al._ [2019] Y.-R. Liu, H.-X. Chen, W. Chen, X. Liu, and S.-L. Zhu, Prog. Part. Nucl. Phys. 107, 237 (2019).
  * Brambilla  _et al._ [2020] N. Brambilla, S. Eidelman, C. Hanhart, A. Nefediev, C.-P. Shen, C. E. Thomas, A. Vairo, and C.-Z. Yuan, Phys. Rep. 873, 1 (2020).
  * Lucha  _et al._ [2021] W. Lucha, D. Melikhov, and H. Sazdjian, Prog. Part. Nucl. Phys. 120, 103867 (2021).
  * Dong  _et al._ [2021] X.-K. Dong, F.-K. Guo, and B.-S. Zou, Commun. Theor. Phys. 73, 125201 (2021).
  * Brambilla [2022] N. Brambilla, PoS LATTICE2021, 020 (2022).
  * Meng  _et al._ [2023] L. Meng, B. Wang, G.-J. Wang, and S.-L. Zhu, Phys. Rep. 1019, 1 (2023).
  * Chen  _et al._ [2023] H.-X. Chen, W. Chen, X. Liu, Y.-R. Liu, and S.-L. Zhu, Rept. Prog. Phys. 86, 026201 (2023).
  * M. Ablikim _et al_. (2023) [BESIII Collaboration] M. Ablikim _et al_. (BESIII Collaboration), Chin. Phys. C 47, 033001 (2023).
  * Aguilar and Combes [1971] J. Aguilar and J. M. Combes, Comm. Math. Phys. 22, 269 (1971).
  * Balslev and Combes [1971] E. Balslev and J. M. Combes, Comm. Math. Phys. 22, 280 (1971).
  * Cheng  _et al._ [2022] J.-B. Cheng, Z.-Y. Lin, and S.-L. Zhu, Phys. Rev. D 106, 016012 (2022).
  * [71] Z.-Y. Lin, J.-B. Cheng, and S.-L. Zhu, arXiv:2205.14628 .
  * P. A. Zyla _et al_. (2020) [Particle Data Group Collaboration] P. A. Zyla _et al_. (Particle Data Group Collaboration), Prog. Theor. Exp. Phys. 2020, 083C01 (2020).
  * Aoyama  _et al._ [2006] S. Aoyama, T. Myo, K. Katō, and K. Ikeda, Prog. Theor. Phys. 116, 1 (2006).
  * Ho [1983] Y. Ho, Phys. Rep. 99, 1 (1983).
  * Ahmed and et al [2001] S. Ahmed and et al, Phys. Rev. Lett. 87, 251801 (2001).
  * Homma  _et al._ [1997] M. Homma, T. Myo, and K. Kato, Prog. Theor. Phys. 97, 561 (1997).
  * Cheng  _et al._ [2023] J.-B. Cheng, D.-X. Zheng, Z.-Y. Lin, and S.-L. Zhu, Phys. Rev. D 107, 054018 (2023).


with the fixed parameters of $f_{\rm{high}-J}\,=\,1$, and $L_{\star}$ = 20 L⊙.
Each column is one radiative transfer run of POLARIS, with
$a_{\rm{max}}\,=\,$0.5, 2, 10, 30, and 50 $\mu$m, while each row corresponds
to a wavelength of observation for the radiative transfer, i.e., 0.87 mm,
1.3mm, and 3mm, as indicated on the first plot of each row. We only show the
polarization fraction.

Figure 13: Comparisons on the evolution of
$\mathcal{S}\times\mathcal{P_{\rm{frac}}}$ as a function of normalized column
density $N_{{\textrm{H}_{2}}}/N_{{\textrm{H}_{2},\textrm{peak}}}$, between
ALMA observations, and our models, before and after spatial filtering. The
left (right) column corresponds to results of the radiative transfer performed
at 0.87 mm (1.3mm). Same as Figure 6 for the sets I and IV (see Table 1).

##  Appendix C The effects of $f_{\textrm{high-J}}$ on the radiative transfer
results

We refer to the recent study by Chau Giang et al. (2022) who thoroughly
explored the polarized dust emission produced by super-paramagnetic grains.
Yet we explore in this appendix the effects of changing the $f_{\rm{high}-J}$
parameter on the results of the radiative transfer calculations. In Figure 14,
we vary the $f_{\rm{high}-J}$ parameter, from 0 to 1, with
$a_{\rm{max}}\,=\,10\,\mu$m and $L_{\star}$ = 20 L⊙. The polarized intensity
and polarization fraction increase with increasing fraction of grains at the
high-$J$ attractor point, the fraction of grains (among those that are
aligned) that can be considered perfectly aligned if their internal alignment
is efficient (Hoang, 2022; Chau Giang et al., 2022). The resulting
polarization fraction is strongly affected by this parameter. Inside the
central $\sim$ 2000 au, $f_{\rm{high}-J}\,=\,0$ causes the polarization
fraction to not be higher than 3%, while $f_{\rm{high}-J}\,=\,1$ causes the
polarization fraction to reach $\sim\,12\%$ in the outflow cavity walls.

Figure 15 presents the evolution of $\mathcal{S}\times\mathcal{P_{\rm{frac}}}$
as a function of normalized column density
$N_{{\textrm{H}_{2}}}/N_{{\textrm{H}_{2},\textrm{peak}}}$, in ALMA
observations and our models, before and after spatial filtering, for the set
II of radiative transfer runs. In this Figure, we vary the
$f_{\textrm{high-J}}$ values in the range 0.0, 0.25, 0.5, 0.75, 1.0., with
$L_{\star}$ = 20 L$\odot$ and $a_{\rm{max}}\,=\,$10 $\mu$m. Only values of
$f_{\textrm{high-J}}$ close to 1 could produce high enough
$\mathcal{S}\times\mathcal{P_{\rm{frac}}}$ values in order to match ALMA
observations.

Figure 14: Effects of the $f_{\rm{high}-J}$ parameter on the radiative
transfer results at 0.87 mm, with the fixed parameters of
$a_{\rm{max}}\,=\,10\,\mu$m, and $L_{\star}$ = 20 L⊙. This correspond to the
set II of Table 1. Each column is one radiative transfer run of POLARIS, with
$f_{\rm{high}-J}\,=\,$0.0, 0.25, 0.5, 0.75, and 1.0. The top row is the
polarized intensity $P$, the bottom row is the polarization fraction.

Figure 15: Comparisons on the evolution of
$\mathcal{S}\times\mathcal{P_{\rm{frac}}}$ as a function of normalized column
density $N_{{\textrm{H}_{2}}}/N_{{\textrm{H}_{2},\textrm{peak}}}$, between
ALMA observations, and our models, before and after spatial filtering. Same as
Figure 6 for the set II (see Table 1).

##  Appendix D Correction for the optical thickness on the normalized column
density of radiative transfer models.

The distribution of normalized gas column density is more peaked in our model
than in the distribution derived from ALMA observations in Le Gouellec et al.
(2020). This difference can be explained by optically thick emission in the
line of sights located toward the total intensity peak in ALMA maps,
precluding us from accessing the true value of the column density peak. This
effect results in the distributions of
$\mathcal{S}\times\mathcal{P_{\rm{frac}}}$ as a function of normalized column
density $N_{{\textrm{H}_{2}}}/N_{{\textrm{H}_{2},\textrm{peak}}}$ for our
models being shifted to the left compared to the ALMA distribution.

We applied a tentative correction to the column density map in our model to
make the comparisons with the ALMA observations more reliable. This correction
factor is computed from the optical depth $\tau$ obtained in each radiative
transfer calculation (see Fig. 2). In each pixel, the gas column density is
multiplied by $(1-e^{-\tau_{\lambda}})/\tau_{\lambda}$, to take into account
the error made in assuming an optically thin emission computing
$N_{{\textrm{H}_{2}}}/N_{{\textrm{H}_{2},\textrm{peak}}}$ in ALMA
observations. This correction is a first order approximation, and represents
the error made assuming an optically thin emission while $\tau$ is no longer
close to 0. Figure 16 presents the radial distribution of $\tau_{\lambda}$,
and the distribution of original and corrected $N_{{\textrm{H}_{2}}}$, at
$\lambda\,=\,1.3$ and 0.87 mm.

Figure 16: Correction for the optical thickness on the gas column density of
our radiative transfer models, computed with $a_{\textrm{max}}$ = 10 $\mu$m,
at $\lambda\,=\,0.87$ mm and $\lambda\,=\,1.3$ mm, respectively. The dot-dash
light (solid dark) blue line shows the azimuthally averaged modified
(original) gas column density as a function of radius. The dark red solid line
shows the azimuthally averaged optical depth $\tau$ as a function of radius.