The properties of the S-wave D_sD_s bound state

Since the discovery of the $ X(3872) $ (also known as $ \chi_{c1}(3872) $) [1], many charmonium-like states have been reported in experiments. However, the masses of these charmonium-like states do not match the excited states predicted by the relativistic quark model [2]. The structures of these charmonium-like states have thus become the focus of interest and provide a unique window for understanding the nature of the strong force. Interestingly, these states share a common feature in that their masses are close to the thresholds of some hadron pairs, which suggests that they may be loose hadronic molecules [3−15].

Not long ago, the LHCb Collaboration reported a new resonant structure, $ X(3960) $, with $ J^{PC}=0^{++} $, observed in the $ D_s^+D_s^- $ mass distribution of the $ B^+\rightarrow D_s^+D_s^-K^+ $ decay with a significance greater than $ 12\sigma $ [16]. The mass and width measured by LHCb are

This structure lies just above the $ D_s^+D_s^- $ threshold and represents an excellent candidate for a $ D_s^+D_s^- $ hadronic molecule. The LHCb Collaboration also reported a scalar state, $ X_0(3930) $, in the $ D^+D^- $ mass distribution of the $ B^+\rightarrow D^+D^-K^+ $ decay [17, 18] in 2020, with a mass and width of

which lies 10 MeV below the $ D_s^+D_s^- $ threshold. The Particle Data Group classifies these as the same particle as the structure first observed by the Belle experiment in the $ \omega J/\psi $ invariant mass spectrum of the $ B \to K \omega J/\psi $ process [19], referred to as X(3915) (or $ \chi_{c0}(3915) $) [20].

There have been extensive theoretical studies of $ X(3915) $, with several investigating its possible assignment as a charmonium state. Within an unquenched quark model, Ref. [25] finds the mass and decay width of $ \chi_{c0}(2P) $ to be consistent with measurements of $ X(3915) $. Using a rescattering mechanism to compute $ {\cal{B}}[B \rightarrow K\chi_{c0}(2P)] $, Ref. [26] argues that the $ \chi_{c0}(2P) $ state plays a crucial role in the LHCb measurement of $ B \rightarrow K D\bar{D} $. Coupled-channel analyses indicate that assigning $ X(3960) $ to the $ \chi_{c0}(2P) $ charmonium state cannot be excluded [27]. Because $ \chi_{c0}(2P) $ predominantly decays to $ D\bar{D} $, Ref. [29] proposes $ e^+e^- \rightarrow \omega D\bar{D} $ as an ideal process to identify this state. In contrast, quark-model calculations employing the $ ^3P_0 $ decay model find that, although the mass of $ \chi_{c0}(2^3P_0) $ is compatible with $ X(3915) $, its predicted strong decay width is much larger than the experimental value [21]. Guo and Meissner [22], analyzing BaBar data, likewise conclude that $ X(3915) $ is difficult to interpret as $ \chi_{c0}(2P) $, a view supported by Ref. [23]. Within a constituent quark framework, Ref. [24] finds $ X(3915) $ to be dominantly molecular, with a bare $ q\bar{q} $ probability below 45%, and favors a $ J^{PC}=2^{++} $ assignment. Finally, an analysis of the $ D\bar{D} $ invariant-mass distribution in $ e^+e^- \rightarrow J/\psi D\bar{D} $ measured by the Belle Collaboration reveals no clear enhancement above the $ D\bar{D} $ threshold to support the existence of a $ \chi_{c0}(2P) $ state [28].

Interpreting X(3915) as a $ D_s\bar{D}_s $ molecular state is another popular approach [8−15, 24, 30−35, 37−40]. By studying the interaction of the $ D\bar{D} $ and $ D_s\bar{D}_s $ coupled channels in $ B^-\rightarrow K^-J/\psi\omega $ decay [8] and $ B^+\rightarrow D^+D^-(D_s^+D_s^-)K^+ $ decay [9], a $ D_s\bar{D}_s $ bound state was found around 3930 MeV, which couples strongly to $ D_s\bar{D}_s $, with no extra resonance signal at 3960 MeV. Such a state was also predicted in lattice QCD simulations [10], with a shallow binding energy of about 6.2 MeV. The authors of Ref. [11] employed an effective field theory based on heavy quark spin symmetry and found that the $ D_s^+D_s^- $ invariant mass distribution data can be well described by either a bound or a virtual state below the $ D_s^+D_s^- $ threshold, with a pole mass of (3928$ \pm $3) MeV. The existence of a $ D_s\bar{D}_s $ bound state from B decay processes was also supported by effective Lagrangian approaches [14]. The $ D_s\bar{D}_s $ system can also form a bound state in a contact-range theory [15]. Based on the vector-meson-dominance model [35], the potential from ϕ exchange is insufficient to form a $ D_s\bar{D}_s $ bound state; instead, a shallow virtual state was obtained with a virtual energy roughly in the range [4.7, 35.5] MeV. No $ D_s\bar{D}_s $ bound state was found in a potential model including ϕ exchange contributions [40]. In Ref. [32], the effects of both ϕ and $ J/\psi $ exchange were considered, and it was found that the $ D_s\bar{D}_s $ system could exist as a bound state within the quasipotential Bethe-Salpeter equation approach. Therefore, elucidating the properties of the $ D_s\bar{D}_s $ bound state is important for understanding the nature of $ X(3930) $ and $ X(3960) $.

The purpose of the present work is to apply the Bethe-Salpeter (BS) equation to investigate the possibility of a $ D_s\bar{D}_s $ bound state when contributions from ϕ and $ J/\psi $ meson exchanges are considered. We will further investigate the partial decay widths of a possible $ D_s\bar{D}_s $ bound state.

This work is organized as follows. In Sec. Ⅱ, we review the basic BS formalism for a system of two pseudoscalar particles. The formalism for the partial decay widths of the $ D_s\bar{D}_s $ bound state to $ D\bar{D} $, $ \eta_c\eta $, and $ J/\psi\omega $ final states will be presented in Sec. Ⅲ. Numerical results are presented in Sec. Ⅳ. A brief summary is given in the final section.

In this section, we derive the BS equation for the S-wave $ D_s\bar{D}_s $ system. The BS wave function for the S-wave $ D_s\bar{D}_s $ bound state with total momentum P is defined as

The coordinate X denotes the center of mass, while p and x represent the relative momentum and relative coordinate of the $ D_s $ and $ \bar{D}_s $ pair, respectively. We define $ \lambda_i = m_{D_s(\bar{D}_s)}/(m_{D_s} + m_{\bar{D}_s}) = 1/2 $ (for $ i = 1,2 $). The individual momenta of $ D_s $ and $ \bar{D}_s $ can then be expressed in terms of P and p as $ p_1 = \lambda_1 P + p $ and $ p_2 = \lambda_2 P - p $, respectively.

The S-wave $ D_s\bar{D}_s $ bound state wave function $ \chi_P(p) $ satisfies the following BS equation:

where $ S_{D_s}(p_1) $ and $ S_{\bar{D}_s}(p_2) $ are the propagators for $ D_s $ and $ \bar{D}_s $, respectively, and $ K_P(P,p,q) $ is the kernel which can be derived from the four-point Green function.

The interaction kernel of the $ D_s\bar{D}_s $ system can be obtained through the exchange of ϕ and $ J/\psi $ mesons. Since the contribution from $ J/\psi $ exchange was found to be significant in the $ D\bar{D}^\ast $ interaction for the production of the $ Z_c(3900) $ [41, 42], and its effect is likewise crucial for the $ D\bar{D} $ system to form a bound state [43], we include the contribution of $ J/\psi $ exchange in the current work. The Lagrangian can be constructed based on chiral symmetry and heavy quark symmetry:

The magnitudes of the coupling parameters are very important for the possibility of forming a bound state. For the coupling constant $ g_{D_sD_s\phi} $, the results from the SU(3) relation ($ g_{D_sD_s\phi}=\dfrac{m_{D_s}}{m_D}g_{DD\rho} $) [44] and from QCDSR [45−47] vary widely, as listed in Table 1. The coupling constant $ g_{D_sD_sJ/\psi} $ has been obtained through flavor SU(3) symmetry ($ g_{D_sD_sJ/\psi} = g_{DDJ/\psi} $) [48], and from the relation with the gauge coupling $ g_2 $ ($ g_{D_sD_sJ/\psi}=2g_2m_D\sqrt{m_{J/\psi}} $ with $ g_2=\sqrt{m_{J/\psi}}/(2m_Df_{J/\psi}) $ and $ f_{J/\psi}=405 $ MeV) [32].

	[44]	[45]	[46]	[47]
$ g_{D_sD_s\phi} $	3.89	$ 1.90^{+0.17}_{-0.13} $	$ 1.45 \pm 0.34 $	$ 2.08^{+0.41}_{-0.43} $

Table 1.  The value of the coupling constant $ g_{D_sD_s\phi} $.

At the tree level, the t-channel interaction kernel for the BS equation in the so-called ladder approximation is

where $ \Delta_{\mu\nu}(p_1-q_1) $ represents the propagator for the ϕ or $ J/\psi $ meson.

Since the strong interaction vertices are determined by the physical particles and the off-shell exchanged particles, it is necessary to introduce a form factor $ F(k) $ to account for the off-shell effects of the t-channel exchanged particles. The form factor is defined as

where k is the momentum of the exchanged particle, $ m_V $ is its physical mass, and Λ is the cutoff in the form factor. The value of the cutoff Λ should be close to the physical mass of the exchanged particle. It can be reparameterized as $ \Lambda=m_V+\alpha \Lambda_{\text{QCD}} $, where $ \Lambda_{\text{QCD}} $ is the QCD scale (approximately 220 MeV) and α is expected to be of order unity. The value of α depends on the exchanged and external particles involved in the strong interaction vertex and cannot be obtained from first principles.

To solve the BS Eq. (4), we use the instantaneous approximation in the kernel. In this approximation, the energy exchanged between the constituent particles of the bound system is neglected. Studies from lattice QCD [10], effective field theory [11], and contact-range theory [15] suggest that the $ D_s\bar{D}_s $ system has a very small binding energy. This indicates that the binding of the constituent particles is weak, making it reasonable to employ the instantaneous approximation in the kernel of the BS equation.

Substituting Eqs. (6) and (7) into Eq. (4) and applying the instantaneous approximation, in the center-of-mass frame of the bound state ($ P=(M,{\bf{0}}) $), we obtain:

for each exchanged meson with mass $ m_V $ and coupling $ g_{D_sD_sV} $, $ \omega_{1(2)} $ is the energy defined as $ \sqrt{m^2_{D_s(\bar{D}_s)}+{\bf{p}}^2} $.

Then, performing the integration over $ p_0 $ on both sides (applying the residue theorem on the right-hand side) and over $ q_0 $ on the right-hand side, Eq. (8) becomes

To solve the integral equation (9), we need to complete the azimuthal integration and discretize the integration region into n pieces (with n sufficiently large). This discretization transforms the integral equation into an eigenvalue equation for the n-dimensional vector $ \chi_P(|{\bf{p}}|) $. Then we can solve the BS equation numerically.

In this section, we apply the obtained BS wave functions to calculate the partial decay widths of the $ D_s\bar{D}_s $ bound state. Prior to this, it is necessary to obtain the normalized BS wave function, which corresponds to determining the physical wave function. The normalization condition for the BS wave function is given by

where $ I_P(p,q)=(2\pi)^4\delta(p-q)S_{D_s}^{-1}(p_1)S_{\bar{D}_s}^{-1}(p_2) $. The energy of the bound state is denoted as $ E_P $, and in the rest frame of the bound state, we have $ P^0=E_P=M $.

After carrying out the $ p_0 $-integrations using appropriate contours and completing the azimuthal integration for the normalization equation (10), we obtain

with

Then, we can obtain the normalized BS wave function. If the $ D_s\bar{D}_s $ bound state wave function obtained in the previous section does not satisfy this normalization equation (the left side of the normalization equation gives some constant $ {\cal{C}}^2\neq 1 $), one must replace $ \chi_P(|{\bf{p}}|)\rightarrow \chi_P(|{\bf{p}}|)/|{\cal{C}}| $ to ensure the correct normalization of the BS wave function.

Now, we proceed to study the decay widths of the $ D_s\bar{D}_s $ bound states. The relevant interaction Lagrangians for the $ D_s\bar{D}_s $ bound state decaying to $ D\bar{D} $, $ \eta_c\eta $, and $ J/\psi\omega $ are [31]

where $ {\cal{D}}=\{{\cal{D}}^{(\ast)0},{\cal{D}}^{(\ast)+},{\cal{D}}_s^{(\ast)+}\} $, and $ {\cal{V}} $ and $ {\cal{P}} $ are the matrices of vector and pseudoscalar mesons, respectively, which have the following forms:

In Eq. (13), the physical ω and ϕ are mixtures of the ideal isospin eigenstates $ \omega^I $ and $ \phi^I $ with the mixing angle $ \theta_v $, where $ \kappa = \cos\theta_v / \sqrt{2} $, $ \zeta = \sin\theta_v / \sqrt{2} $, $ \lambda = -\sin\theta_v $, and $ \sigma = \cos\theta_v $. Similarly, in Eq. (14), η and $ \eta' $ are mixtures of the ideal isospin eigenstates $ \eta^I $ and $ \eta^{'I} $ with the mixing angle $ \theta_p $, where $ \tau = (\cos\theta_p - \sqrt{2}\sin\theta_p) / \sqrt{6} $, $ \xi = (\sin\theta_p + \sqrt{2}\cos\theta_p) / \sqrt{6} $, $ \gamma = (-2\cos\theta_p - \sqrt{2}\sin\theta_p) / \sqrt{6} $, and $ \delta = (-2\sin\theta_p + \sqrt{2}\cos\theta_p) / \sqrt{6} $. In this work, we adopt the central value $ \theta_v = (3.4 \pm 0.2)^\circ $, as determined in Refs. [49−51]. We also take $ \theta_p = -19.1^\circ $ as in Refs. [52, 53]. The relevant coupling constants can be estimated using the heavy quark limit and chiral symmetry.

The relevant coupling constants in Eq. (15) of our manuscript can be estimated in the heavy quark limit and with chiral symmetry, where the coupling constants of the S-wave charmonia and charmed mesons can be determined by

where the gauge coupling $ g_2 $ relates to the decay constant of $ J/\psi $ by $ g_2=\sqrt{m_{J/\psi}}/(2m_Df_{J/\psi}) $, which is estimated from the leptonic decay width of $ J/\psi $ [31]. The other coupling constants can be estimated by

where $ \beta_V $ can be fixed to 0.9 by vector meson dominance, $ \lambda_V $ is obtained from light-cone sum rules and lattice QCD as $ 0.56\, {\rm{GeV}}^{-1} $. Imposing the Kawarabayashi-Suzuki-Riazuddin-Fayyazuddin relations, $ g_V = m_\rho/f_\pi $, with $ f_\pi = 0.132 $ GeV (the pion decay constant), and $ g = 0.59 $ (estimated from the partial width for $ D^\ast\rightarrow D\pi $) [56]. Since we have accounted for variations in the coupling constant $ g_{D_sD_s\phi} $ from different approaches while solving the BS equation—and these variations span a wide range—we proceed with the analysis.

With the effective Lagrangians in Eq. (12), the lowest-order amplitudes for the decays of the $ D_s\bar{D}_s $ bound state to $ D\bar{D} $, $ \eta_c\eta $, and $ J/\psi\omega $ are

respectively, where $ p'_i (i=1,2) $ is the momentum of the i-th particle in the final state. In the rest frame of the bound state, the momenta of the two particles in the final state can be taken as $ p'_1=(E_1,{\bf{p}}') $, $ p'_2=(E_2,-{\bf{p}}') $ with $ E_1=\sqrt{{\bf{p}}^{'2}+m_1^2} $ and $ E_2=\sqrt{{\bf{p}}^{'2}+m_2^2} $. We define $ p'\equiv \lambda_2 p'_1-\lambda_1p'_2 $, which is not the relative momentum of the final particles; then, $ p'=(\lambda_2E_1-\lambda_1E_2,{\bf{p}}') $.

In the rest frame, the two-body decay width of the bound state can be written as

where

the Källén function is given by $ \lambda(a,b,c)=a^2+b^2+c^2- 2ab-2ac-2bc $.

The coupling constants and form factors associated with strong interaction vertices play a crucial role in studying possible bound states of the $ D_s\bar{D}_s $ system and in estimating partial decay widths of the $ D_s\bar{D}_s $ bound state into $ D\bar{D} $, $ \eta_c\eta $, and $ J/\psi\omega $ final states. As shown in Table 1, the magnitude of $ g_{D_sD_s\phi} $ varies considerably depending on whether it is obtained from the SU(3) relation or from QCD sum rules. In our work, we consider the range of $ g_{D_sD_s\phi} $ values from the minimum to the maximum presented in Table 1.

Since the exchanged particle is not on-shell, it is necessary to introduce a form factor $ F(k) $ to account for the off-shell effect of the t-channel exchanged particle. The cutoff Λ in the form factor $ F(k) $ should be close to the physical mass of the exchanged particle. We take $ \Lambda=m+\alpha\Lambda_{\text{QCD}} $, where α depends on the exchanged particle and the external particles involved in the strong interaction vertex. The parameter α is expected to be of order unity and cannot be determined from first principles. In this work, we identify all possible values of α for which the $ D_s\bar{D}_s $ system can form a bound state.

We then solve the BS equation (9) numerically by discretizing it into a matrix eigenvalue equation using the Gaussian quadrature method. For each trial set of values of the coupling constant $ g_{D_sD_s\phi} $, the parameter α, and the binding energy $ E_b $ of the $ D_s\bar{D}_s $ system (defined as $ E_b=m_{D_s}+m_{\bar{D}_s}-E $ and varying from 0 to 15 MeV), we obtain all eigenvalues of this equation. The set of values of $ g_{D_sD_s\phi} $, α, and $ E_b $ that yields the eigenvalue closest to 1.0 is selected. We find that for a very small coupling constant $ g_{D_sD_s\phi} $, a very large parameter α is required for a bound state to exist. Therefore, we choose $ g_{D_sD_s\phi} $ in the interval of (2–3.89) for our analysis. The values of $ g_{D_sD_s\phi} $ and α for possible bound states with binding energies $ E_b=1 $ MeV, $ E_b=5 $ MeV, and $ E_b=15 $ MeV are presented in Fig. 1. From Fig. 1, we observe that for a given coupling constant $ g_{D_sD_s\phi} $, a larger binding energy requires a larger value of the parameter α. Additionally, our study reveals that the $ D_s\bar{D}_s $ system cannot form a bound state if only the contribution of ϕ or $ J/\psi $ exchange is considered.

Figure 1.  (color online) The possible bound states in the $ D_s\bar{D}_s $ system.

Based on the obtained numerical results for the $ D_s\bar{D}_s $ bound state wave functions and the corresponding parameter α, we calculate the partial decay widths of the $ D_s\bar{D}_s $ bound state. The uncertainties from the coupling constants involved in the decay processes are considered. These coupling constants are derived via SU(3)-flavor symmetry, with symmetry breaking at the level of $ \delta = 19 $% [14]. Following the approach in Ref. [54], we estimate the uncertainty in the partial decay widths using $ \Gamma = \Gamma(1+\delta)^2 $. The partial decay widths for the $ D_s\bar{D}_s $ bound state decaying into $ D\bar{D} $, $ \eta_c\eta $, and $ J/\psi\omega $ final states are shown in Figs. 2, 3, and 4, respectively. Our results show that these decay widths vary with the parameter α and are sensitive to α in some regions. Moreover, for a fixed value of α, the partial decay widths of the $ D_s\bar{D}_s $ bound state increase as the binding energy increases.

Figure 2.  (color online) The particle decay width of the $ D_s\bar{D}_s $ bound state to $ D\bar{D} $ is shown for $ E_b $ = 1 MeV (a), $ E_b $ = 5 MeV (b), and $ E_b $ = 15 MeV (c) across the allowed parameter α.

Figure 3.  (color online) The particle decay width of the $ D_s\bar{D}_s $ bound state to $ \eta_c\eta $ is shown for $ E_b $ = 1 MeV (a), $ E_b $ = 5 MeV (b), and $ E_b $ = 15 MeV (c) across the allowed parameter range of α.

Figure 4.  (color online) The decay width of the $ D_s\bar{D}_s $ bound state to $ J/\psi\omega $ is shown for $ E_b $ = 1 MeV (a), $ E_b $ = 5 MeV (b), and $ E_b $ = 15 MeV (c), across the allowed range of the parameter α.

From Figs. 2, 3, and 4, we observe that the $ D\bar{D} $ channel dominates the decays of the $ D_s \bar{D}_s $ molecule, while the $ J/\psi\omega $ width is the smallest. This hierarchy occurs because $ D\bar{D} $ proceeds via light $ K^\ast $ exchange, whose propagator is less suppressed off shell; in contrast, the $ \eta_c \eta $ and $ J/\psi\omega $ modes require heavy $ D_s^{(\ast)} $ exchange, which introduces heavier propagators and stronger suppression. We also include η-$ \eta' $ and ω-ϕ mixing, with mixing angles $ \theta_p = -19.1^\circ $ and $ \theta_v = (3.4\pm0.2)^\circ $, which further suppress these two channels. The $ J/\psi\omega $ channel is particularly suppressed due to the tiny $ s\bar{s} $ component in ω and the small $ g_{D_s^{(\ast)}D_s\omega} $ coupling constants, rendering its width orders of magnitude smaller. Additionally, for very shallow binding (1–5 MeV), we observe that $ \Gamma(\eta_c \eta) $ does not increase monotonically with α. This may be because the two Feynman diagrams ($ \eta_c $ emitted from the $ D_s $ versus the $ \bar{D}_s $ vertex) exhibit different α dependences. However, it is difficult to attribute this behavior purely to kinematics or dynamics. Importantly, this nonmonotonicity occurs only for extremely shallow binding in a subleading channel and does not affect the overall decay pattern or our main conclusions.

In Ref. [31], the authors also investigated the decay behaviors of the $ D_s^+D_s^- $ bound state using an effective Lagrangian approach. Our results for the decay widths to $ D\bar{D} $ and $ J/\psi\omega $ final states are consistent with those obtained in Ref. [31]. However, there is a discrepancy in the decay width to the $ \eta_c\eta $ final state. In Ref. [31], the decay width for $ \eta_c\eta $ ranges from over 1 MeV to nearly 100 MeV as the parameter Λ varies from 0.2 to 1 GeV. In contrast, in our study, the decay width for $ \eta_c\eta $ varies by only a few MeV within the allowed range of the parameter α.

Furthermore, the $ X(3915) $ state, initially discovered by the Belle Collaboration [19] and subsequently assigned $ J^{PC} = 0^{++} $ by the BaBar [55] and LHCb [18] Collaborations, has been considered a possible candidate for a charmonium-like state [22, 23]. The mass of $ X(3915) $ lies just below the $ D_s\bar{D}_s $ threshold by about 16 MeV [20]. Ref. [30] suggested that the properties of $ X(3915) $ can be explained if it is an S-wave $ D_s\bar{D}_s $ bound state with a binding energy of about 18 MeV. In our model, using the average mass of $ X(3915) $ from the Particle Data Group [20], 3921.7 MeV, we find the partial decay widths to the $ D\bar{D} $, $ \eta_c\eta $, and $ J/\psi\omega $ final states to be in the ranges (35.06–135.86) MeV, (0.51–3.04) MeV, and (0.49–2.05)$ \times10^{-3} $ MeV, respectively. For the shallow $ D_s\bar{D}_s $ bound state with $ E_b = 6.2 $ MeV predicted by lattice QCD, the partial decay widths to $ D\bar{D} $, $ \eta_c\eta $, and $ J/\psi\omega $ final states are (21.92–94.77) MeV, (0.37–1.40) MeV, and (0.05–0.69)$ \times10^{-3} $ MeV, respectively.

In this work, we have investigated possible bound states of the $ D_s\bar{D}_s $ system using the BS equation and calculated the partial decay widths of these states to the $ D\bar{D} $, $ \eta_c\eta $, and $ J/\psi\omega $ final states. To accomplish this, we applied the ladder and instantaneous approximations to the kernel, which includes contributions from ϕ and $ J/\psi $ exchange diagrams. Our calculations show that the $ D_s\bar{D}_s $ system can form bound states when both ϕ and $ J/\psi $ exchange contributions are considered. However, when only one of these exchanges is included, no bound states are formed. The possible bound states of the $ D_s\bar{D}_s $ system for binding energies $ E_b $ = 1, 5, and 15 MeV are presented in Fig. 1.

We used the numerical solutions of the BS wave functions to compute the partial decay widths of the $ D_s\bar{D}_s $ bound states into $ D\bar{D} $, $ \eta_c\eta $, and $ J/\psi\omega $ final states. The results for these partial decay widths are shown in Figs. 2, 3, and 4, respectively. Our results indicate that the partial decay widths depend on the parameter α and are sensitive to its value in certain regions. We observe that the dominant decay channel for the $ D_s\bar{D}_s $ bound state is $ D\bar{D} $. Therefore, neglecting contributions from other channels such as light meson decays and photon radiation, the total decay width of the $ D_s\bar{D}_s $ bound state is estimated to be in the ranges 8.56–46.39 MeV, 20.12–102.43 MeV, and 35.83–174.12 MeV for binding energies $ E_b $ = 1, 5, and 15 MeV, respectively. These results show that the decay width increases with binding energy. Based on our calculations, the values of α, and the magnitude of the decay widths, it is plausible to conclude that the $ D_s\bar{D}_s $ system can form a loosely bound state. We hope that our theoretical predictions regarding the properties of the $ D_s\bar{D}_s $ bound state will be tested in future experimental studies.

One of the authors (Zhen-Yang Wang) would like to thank Professors Feng-Kun Guo and Bing-Song Zou for their financial support during his visit to the Institute of Theoretical Physics, Chinese Academy of Sciences, as well as for helpful discussions with Dr. Pan-Pan Shi.