$DK$, $DDK$, and $DDDK$ molecules--understanding the nature of the $D_{s0}^*(2317)$
Tian-Wei Wu Ming-Zhu Liu Li-Sheng Geng Emiko Hiyama Manuel Pavon Valderrama



School of Physics and Nuclear Energy Engineering,
  Beihang University, Beijing 100191, China


School of Physics and Nuclear Energy Engineering,
  Beihang University, Beijing 100191, China

    
    lisheng.geng@buaa.edu.cn
    School of Physics and Nuclear Energy Engineering, Beihang University, Beijing 100191, China
    Beijing Key Laboratory of Advanced Nuclear Materials and
      Physics, Beihang University, Beijing 100191, China
      School of Physics and Engineering, Zhengzhou University, Zhengzhou, Henan 450001, China
      
    
    hiyama@riken.jp
    Department of Physics, Kyushu University, Fukuoka 819-0395,
      Japan, RIKEN Nishina Center, RIKEN, Wako 351-0198, Japan
      Nishina Center for Accelerator-Based Science, RIKEN, Wako, 351-0198, Japan



    
    mpavon@buaa.edu.cn
    School of Physics and Nuclear Energy Engineering,
      Beihang University, Beijing 100191, China



The DK interaction is strong enough to form a bound state,
the D_s0^*(2317).

This in turn begs the question of whether there are bound states
composed of several charmed mesons and a kaon.

Previous calculations indicate that the three-body DDK system
is probably bound, where the quantum numbers are J^P = 0^-,
I=12, S = 1 and C = 2.

The minimum quark content of this state is ccq̅s̅ with q=u,d,
which means that, if discovered, it will be an explicitly exotic tetraquark.

In the present work. we apply the Gaussian Expansion Method to study the DDDK system and show that it binds as well. 

The existence of these three and four body states is rather robust
with respect to the DD interaction and subleading (chiral)
corrections to the DK interaction.


If these states exist, it is quite likely that their heavy quark symmetry counterparts exist as well. 

These three-body DDK and four-body DDDK molecular states could be viewed as counterparts of  atomic nuclei, which are clusters of nucleons bound by the residual strong force, or chemical molecules, which are clusters of atoms bound by the residual electromagnetic interaction. 





DK, DDK, and DDDK molecules–understanding the nature of
the D_s0^*(2317)
    Manuel Pavon Valderrama
    December 30, 2023
========================================================================






§ INTRODUCTION

In 2003 the BaBar collaboration discovered
the D_s0^*(2317) <cit.> [From now on, we will
simply refer to it as D_s0^* unless specified otherwise.], a strange-charmed scalar meson,
the observation of which was subsequently confirmed
by CLEO <cit.> and Belle <cit.>.

Its mass is about 160  MeV below the one predicted
for the lightest cs̅ scalar state in the naive quark model,
which makes it difficult to interpret the D_s0^*  
as a conventional qq̅ state <cit.>.

On the other hand, the D_s0^* can be easily explained as a
dynamically generated state arising from the Weinberg-Tomozawa (WT)
DK interaction <cit.>.

This has led to the prevailing idea that the D_s0^*(2317)
is a molecular state, a hypothesis which has been further supported by
a series of Lattice QCD
simulations <cit.>. For a recent brief summary of all the experimental, lattice QCD, and theoretical supports for such
an assignment, see, e.g., Ref. <cit.>.

If the DK interaction is strongly attractive, a natural question to ask is
what happens when one adds one extra D meson to
the system [It has been shown that the DD̅^*K system
  binds as well in two recent works <cit.>, though
  the dynamics in these two frameworks are quite different.].

The answer seems to be that it binds <cit.>.

In Ref. <cit.> it was noticed that the D D_s0^*
system can exchange a kaon near the mass shell, leading to a relatively
long-range attractive Yukawa potential that is strong enough to bind.

This conclusion is left unchanged if one explicitly considers the composite
nature of the D_s0^*, which simply leads
to more binding <cit.>.

A later, more complete calculation in Ref. <cit.>
leads to a binding energy of about 90 MeV
for the DDK three-body system.


In the present manuscript we revisit the calculation of the DDK bound state
and extend it to the DDDK system by using the Gaussian Expansion Method (GEM),
which  offers a number of advantages compared to previous
studies <cit.>.

First, it allows one to calculate directly the density distribution of
the three (four) body system, which then gives a transparent picture
for their spacial distributions.

Second, it has enough flexibility so that one can study the impact of
the existence of a repulsive core.

Indeed, the chiral potential kernel up to the next to leading order
with the low-energy constants determined by the corresponding
lattice QCD data shows that this may indeed
be the case <cit.>.

The outcome of the exploration presented in this work is that
both the DDK and DDDK systems bind, with binding energies of
the order of 65-70 and 90-100  MeV in each case.

While the DDK bound state, owing to its cc q̅s̅ quark content,
might be produced in experiments in the future, the DDDK bound state
is more likely to be observed on the lattice instead.

This article is organized as follows.

In Sec. <ref>, we explain how we parametrize and determine
the two-body DK and DD interactions.

In Sec. <ref>, we explain how to construct the three- and four-body
DDK and DDDK wave functions and solve the corresponding
Schrödinger equation using the GEM.

In Sec. <ref>, we present our predictions
for the DDK and DDDK bound states and discuss their sensitivity
to a series of possible corrections.

Finally, we summarize the results of this manuscript
in Sect. <ref>.



§ THE S-WAVE DK AND DD POTENTIALS





The calculation of the DDK and DDDK bound states depends on the DK and DD
two-body interactions.

While the DK interaction can be well constrained directly from the assumption
that the D_s0^*(2317) is a DK bound state, and indirectly
from chiral perturbation theory, the DD interaction
is far from being well determined and we will have
to resort to phenomenological models instead.

In this section we will explain the type of potentials we will use to model
these two-body interactions.




 §.§ The DK interaction


The most important contribution to the DK interaction is
the WT term between a D meson and a kaon [Coupled channel
interactions are small, see, e.g. Ref. <cit.>. Therefore, we would work in the single-channel scenario.].

In the non-relativistic limit we can write this interaction
as a standard quantum mechanical potential,

    V_DK(q⃗) = -C_W(I)/2 f_π^2

where the pion decay constant f_π≈ 130 MeV and C_W(I) represents the strength of the WT interaction,
which is

    C_W(0) = 2      C_W(1) = 0   ,

depending on whether we are considering the isospin I=0 or I=1
configuration of the DK system.

The Fourier-transform of the previous potential in coordinate space is

    V_DK(r⃗) = -C_W(I)/2 f_π^2 δ^(3)(r⃗)   ,

which has to be regularized before being used within the Schrödinger equation.

A possible choice is to use a local Gaussian regulator of the type

    V_DK(r; R_c) = -C_W(I)/2 f_π^2 e^-(r/R_c)^2/π^3/2 R_c^3  ,

where R_c is the cutoff we use to smear the delta function.

For sensible choices of the cutoff, this potential reproduces
the D_s0^* pole.

Nowadays we consider the WT interaction
as the leading order (LO) term in the chiral expansion of
the DK potential <cit.>.

In this regard it is interesting to notice that even though
LO chiral perturbation theory (ChPT) indeed indicates that
the I = 0 DK interaction in S-wave is attractive,
it happens that the next-to-leading order (NLO) correction is weakly repulsive,
see e.g. Ref. <cit.>.

This motivates the inclusion of a short-range repulsive core
in the DK interaction, as we will explain in the next paragraph.

For the present purposes a more practical approach will be to consider
the DK interaction in a contact-range effective field theory,
in which at LO we have the (already regularized) potential

    V_DK(r; R_c) = C(R_C) e^-(r/R_c)^2/π^3/2 R_c^3  ,

with R_c the cutoff and where the C(R_c) is now a running coupling constant.

The differences with a unitarized WT term are
(i) that we let the cutoff R_c to float and
(ii) that we consider the strength of the interaction to run with the cutoff.

In this way by varying the cutoff within a sensible range, for which we choose
R_c = 1-3  fm in this work, we can estimate the uncertainty
in the calculations coming from subleading corrections.

We advance that the cutoff variation will be tiny.

Besides the variation of the cutoff,
we will consider a second method to assess the error in our calculations.

Inspired by the fact that ChPT predicts a repulsive core in the DK interaction
at NLO (as previously mentioned), we can explicitly include
this core in the potential

    V_DK(r⃗; R_c) = C_S  e^-(r/R_S)^2/π^3/2 R_S^3 +
      C(R_C) e^-(r/R_c)^2/π^3/2 R_c^3=C'_S e^-(r/R_S)^2+C'_L e^-(r/R_c)^2  ,

where C_S is a coupling constant that we set as to provide a repulsive core,
i.e. we take C'_S> |C'_L|, and R_S is a second cutoff which fulfills
the condition R_S < R_c. For concreteness we take R_S = 0.5  fm.



 §.§ The DD interaction



The DD interaction is not known experimentally,
but there are phenomenological models for it.

Here we will consider the one boson exchange (OBE) potential, which provides
a very simple and intuitive description of the hadron-hadron interactions.

The first qualitatively successful description of the two-nucleon potential
used the OBE model <cit.>,
and the same is true for the first speculations
about the existence of heavy hadron
molecules <cit.>.

The particular version of the OBE model that we will use is the one
in Ref. <cit.>, developed for the description of
heavy meson-meson and heavy meson-antimeson systems.

In the particular case of the DD two-body system, the OBE potential involves
the exchange of the σ, ρ and ω mesons:

    V_DD(r; Λ)=
      V_ρ(r; Λ) + V_ω(r; Λ) + V_σ(r; Λ)

where the contribution of each light meson is regularized by means of
a form factor and Λ is a cutoff.

The particular contribution of each meson can be written as <cit.>

    V_σ(r; Λ)    =   
      - g_σ^2   m_σ 
      W_C(m_σ r, Λ/m_σ )   , 
     
      V_ρ(r; Λ)    =   
      + τ⃗_1 ·τ⃗_2   g_ρ^2   m_ρ 
      W_C(m_ρ r, Λ/m_ρ )   , 
    
      V_ω(r; Λ)    =   
      + g_ω^2   m_ω 
      W_C(m_ω r, Λ/m_ω )   ,

where

    W_C(x, λ) = e^-x/4π x -
      λ e^-λ x/4πλ x -
      (λ^2 - 1)/2λ e^-λ x/4π.

The masses of the bosons we use are m_ρ=0.770 GeV, m_ω=0.780 GeV,
m_σ=0.6 GeV, and the couplings are g_ρ = g_ω = 2.6,
g_σ = 3.4.

The cutoff is set by reproducing the X(3872) pole,
yielding Λ = 1.01^+0.19_-0.10  GeV <cit.>.

Here for the sake of simplicity we will set the cutoff to Λ = 1.0  GeV,
where we note that the cutoff dependence is weak.











§ GAUSSIAN EXPANSION METHOD TO SOLVE THE 3-BODY DDK
   AND 4-BODY DDDK SYSTEMS 



In this section we briefly explain the
Gaussian Expansion Method (GEM) <cit.>
as applied to the DDK and DDDK systems.

In the past the GEM has been successfully applied in hypernuclear
as well as heavy-hadron systems.

The focus of the manuscript is on the one hand to confirm the previous
theoretical studies about the existence of a DDK bound state and
to explore whether there are also bound DDDK tetramers.

Regarding the DDK system, it was investigated
in Ref. <cit.> first
as a D D_s0^* two-body system, a description which is valid
provided that the size of the DDK trimer is larger than its components
(in particular the D_s0^* meson), and second as a genuine
three-body system by solving the Faddeev equations.

In each case the bound state is at about (50-60)  MeV

and (60-100)  MeV 
below the DDK threshold, respectively.

Later a more complete study appeared in Ref. <cit.>,
which uses the method developed by the Valencia group <cit.> to solve the Faddeev equation <cit.>
for the DDK system, predicting a bound state
at about 90  MeV below the DDK threshold.




 §.§ Three-body DDK system


The Schrödinger equation of the DDK 3-body system is
 
    HΨ_JM^total=EΨ_JM^total,

with the corresponding Hamiltonian

    Ĥ=∑_i=1^3p_i^2/2m_i-T_c.m.+∑_1=i<j^3V(r_ij),

where T_c.m. is the kinetic energy of the center of mass and
V(r_ij) is the potential between the i-th and
the j-th particle pair.

The three Jacobi coordinates for the DDK system are shown in Fig. <ref>.

The total wave function is a sum of the amplitudes of the three possible
rearrangement of the Jacobi coordinates, i.e. of the channels (c = 1-3)
shown in Fig. <ref>
  
    Ψ_JM^total=∑_c,αC_c,α Ψ_JM,α^c(𝐫_c,𝐑_c)   ,

  where α={nl,NL,Λ,tT} and C_c,α are the expansion
  coefficients.
  
  Here l and L are the orbital angular momenta for the coordinates
  r and R,  t is the isospin of the  two-body subsystem in each channel,
  Λ and T are the total orbital angular momentum and isospin,
  n and N are the numbers of Gaussian basis function corresponding
  to coordinates r and R, respectively.
  
  For the DD and DK two-body potentials we refer to Sect. <ref>.
  
  The eigen energy E and coefficients are determined by
  the Rayleigh-Ritz variational principle.
  
  Considering that the two D mesons are identical, the total wave function
  should be symmetric with respect to the exchange of the two D mesons,
  which requires that

    P_12Ψ_JM^total=  Ψ_JM^total,

and P_12 is the exchange operator of particles 1 and 2.
The wave function of  each channel has the following form
  
    Ψ_JM,α^c(𝐫_c,𝐑_c) =
        H_T,t^c⊗[Φ_lL,Λ^c]_JM  ,

where H_T,t^c is the isospin wave function, and Φ_lL,Λ^c the spacial wave function. The total isospin wave function reads as
 
 
 
      
    H_T,t^c=1    =[[η_1/2(D_2)η_1/2(K_3)]_t_1η_1/2(D_1)]_1/2  , 
    
              H_T,t^c=2    =[[η_1/2(D_1)η_1/2(K_3)]_t_2η_1/2(D_2)]_1/2  , 
    
              H_T,t^c=3    =[[η_1/2(D_1)η_1/2(D_2)]_t_3η_1/2(K_3)]_1/2  ,

      where η is the isospin wave function of each particle.
      
      The spacial wave function Φ_lL,Λ^c is given in terms of
      the Gaussian basis functions
      
    Φ_lL,Λ^c(𝐫_c,𝐑_c)=[ϕ_n_cl_c^G(𝐫_c)ψ_N_cL_c^G(𝐑_c)]_Λ,

      
    ϕ_nlm^G(𝐫_c)=N_nlr_c^le^-ν_n r_c^2 Y_lm(r̂_c)   ,

  
    ψ_NLM^G(𝐑_c)=N_NLR_c^Le^-λ_n R_c^2 Y_LM(R̂_c)   .

  Here N_nl(N_NL) are the normalization constants of the Gaussian basis
  and the range parameters ν_n and λ_n are given by
  
    ν_n    =1/r_n^2,     r_n=r_mina^n-1   (n=1,n_max)   , 
    λ_N    =1/R_N^2,   R_N=R_minA^N-1   (N=1,N_max)   ,

  in which {n_max,r_min,a or r_max} and  {N_max,R_min,A or
  R_max} are Gaussian basis parameters.
  
  After the basis expansion, the Schrödinger equation of this system
  is transformed into a generalized matrix eigenvalue problem:

    [T_αα'^ab+V_αα'^ab-EN_αα'^ab] 
      C_b,α' = 0
        .

Here, T_αα'^ab is the kinetic matrix element, V_αα'^ab is the potential matrix element and N_αα'^ab is the normalization matrix element.

The quantum numbers of all the allowed configurations are determined by angular momentum conservation, isospin conservation, parity conservation, and Bose-Einstein statistics.  Given that
we only consider S-wave interactions, and only the DK interaction in I=0 is dominant, we obtain the allowed configurations shown in Table <ref>. The DDK system that we are interested in has isospin 1/2 and spin parity 0^-. 


















































 §.§ Four-body DDDK system


A generic four-body system has 18 Jacobi coordinates.

In the DDDK system, owing to the fact that there are three identical
D mesons, the possible configurations of the Jacobi coordinates
reduce to three K-type channels and one H-type channel,
see Fig.[<ref>].

There are 4 identical Jacobi coordinates for each K-type channel
and 6 identical Jacobi coordinates for the H-type channel.

The total wave function of this DDDK system is 
  
    Ψ_I(J^P)^total=∑_c,αA_c,αΨ_α^c(r_c,R_c,ρ_c),      c=1-18   ,

and the wave function in each Jacobi channel reads
  
    Ψ_α^c(r_c,R_c,ρ_c)=H_t,T,I^c⊗Φ_lLλ,σΛ^c,JP  .

  Here t,T,I are the isospin of the coordinates r, R and ρ
  in each channel; l, L and λ are the orbital angular momenta
  for the coordinates r,R and ρ, while σ is the coupling
  of l and L, Λ is the coupling of σ and λ,
  and J, P is the total angular momentum and parity.
  
  The Gaussian basis and parameters are in the same form as those
  in the  3-body system, which are
    
    Φ_lLλ,σΛ^c=[ϕ_n_cl_c^G(r_c)ψ_N_cL_c^G(R_c)]_σ_cφ_ν_cλ_c^G(ρ_c)]_Λ  ,


    ϕ_nlm^G(r_c)=N_nlr_c^le^-ν_n r_c^2 Y_lm(r̂_c)   ,


    ψ_NLM^G(R_c)=N_NLR_c^Le^-λ_N R_c^2 Y_LM(R̂_c)   ,


    φ_νλμ^G(ρ_c)=
      N_νλρ_c^λe^-ω_νρ_c^2 Y_λμ(ρ̂_c)   .

Here N_nl(N_NL) are the normalization constants of
the Gaussian basis and the range parameters
ν_n, λ_n and ω_ν
are given by

    ν_n    =1/r_n^2,     r_n=r_mina^n-1   (n=1,n_max)   , 
    λ_N    =1/R_N^2,   R_N=R_minA^N-1   (N=1,N_max)   , 
    ω_ν   =1/ρ_ν^2,    ρ_ν=ρ_minα^ν-1   (ν=1,ν_max)   .

Since we are considering only S-wave interactions,
we have J=l=L_λ=σ=Λ=0, and the parity is +.

The procedure to  determine the allowed configurations for the DDDK system
is the same as the DDK case.

The  4-body DDDK configurations are shown in Table.<ref>.




§ PREDICTIONS



In this section we discuss the predictions
we make for the DDK and DDDK bound states.

With the two-body inputs of Sect. <ref> and
the three(four)-body configurations detailed in Sect. <ref>,
we can predict the existence of DDK and DDDK bound states.

The outcome is that the DDK trimer will bind by about 70  MeV
and the DDDK tetramer by about 100  MeV, with variations
of a few MeV at most, stemming from the uncertainties
in the DK and DD potentials.



 §.§ Solving the DDK and DDDK systems


The two basic input blocks for the calculation of the DDK and DDDK systems
are the DK and DD interactions, of which the DK one is
the most important factor when it comes to binding.

The DK potential contains the running coupling C(R_c) and the cutoff R_c,
where R_c = (1-3)  fm and C(R_c) is determined from the condition
of reproducing the well-known D_s0^*(2317) as a DK bound state
with a binding energy of 45  MeV.

In addition there are two additional parameters, the coupling C_S
and the short-range radius R_s = 0.5  fm,
which are used to estimate the uncertainties
in the DK potential.

We study three combinations of R_S and R_c, which can be consulted
in Table <ref>, where we also list the values of
the couplings C_R and C(R_c) and the binding energies
of the DDK and DDDK systems.

The different potentials investigated are shown in Fig. <ref>
and the probability density distributions of the DK pair corresponding to
the potentials are shown in Fig. <ref>.








A few comments about the results of Table <ref> are in order.

The first thing we notice is that the impact of the DD interaction is mild.

It makes the DDK and DDDK systems more bound, but only
by a few MeV.

This is a bit relieving as the DD interaction is not well known.

The second interesting observation is that the existence of
the DDK and DDDK bound states is rather robust
with respect to the likely existence of
a short-range repulsive core.

In other words, the existence of the DDK and DDDK bound states is
almost guaranteed as long as the D_s0^* is dominantly
a DK bound state (we will later check that this will still
be the case even if the D_s0^* is a compact
c̅ s state).

The third observation is that as the range of the attraction becomes larger,
two bound state solutions appear instead of one, with the deepest
bound one becoming slightly shallower. 





In Table <ref> we show the root mean square (RMS) radius of
the DK and DDK systems as well as the expectation values
of the kinetic and potential terms.

The RMS radius of the D_s0^*, which ranges from 1.2 to
2.6 fm, increases with the cutoff R_c and
with the coupling C_S of the short-range
repulsive core.

In the DDK system, the RMS radius of the DK pair is slightly
larger than its counterpart in the D_s0^*.

The RMS radius of the DD system also increases if we increase the cutoff
R_c or the coupling C_S.

We notice that the geometry of the DDK system is more or less of
a proper triangle, which agrees qualitatively
with the findings of Ref. <cit.>.

From the last two columns of Table <ref>,
it is clear that the DD interaction is weakly attractive,
accounting for only a few MeV of the total potential energy.
























 §.§ Solving the DDDK system as an equivalent  DDD_s0^* system 


If the separation of the DK pair within the DDK trimer and DDDK tetramer
is comparable to or larger than the expected size of the D_s0^*,
in a first approximation it will be possible to treat the D_s0^*
as a point-like particle, with its compound structure providing
subleading corrections to this point-like approximation.

From Table <ref> we can see that the RMS of the DK subsystem
in the DDK and DDDK systems is similar to that of
the D_s0^* as a DK molecule.

In this regard we notice that in Ref. <cit.>
the D_s0^* is approximated as point-like, where the interaction
between the D and D_s0^* is mediated by one kaon exchange
and is strong enough to form a bound state.

This D D_s0^* molecule is predicted to be 50-60  MeV below
the DDK threshold, to be compared with 65  MeV when we
consider it as a genuine DDK three-body state and ignore
the DD interaction (see Table <ref>).

This indicates that the predictions of the point-like approximation are
reasonably good (for such a simple approximation) and that
the compound structure of the D_s0^* provides
additional attraction.

In the following lines we will extend the ideas of Ref. <cit.>
to the DDDK tetramer, i.e. we will treat it as a three-body DDD_s0^*
system where the D_s0^* is assumed to be a compact meson.

To do this, we first reproduce the two-body calculation
of Ref. <cit.>, but in coordinate space,
and then study the three-body DDD_s0^* system
using the GEM.

The interaction of DD_s0^* is attractive and
reads as
    V_OKE(q⃗)=-h^2ω_K^2/f_π^21/μ_K^2+q⃗^2,

where ω_K=m_D_s0^* -m_D and the effective kaon mass
μ_K=√(m_K^2-ω_K^2).

As in Ref. <cit.>, we take h=0.7 and f_π=130 MeV.

We regularize the potential by multiplying it with a dipole form factor of
the type:

    F_D(q^2) =(Λ^2-m_K^2)^2/(Λ^2-q^2)^2  .

After the inclusion of this form factor,
the DD_s0^* potential in coordinate space reads

    V_DD_so^*(r)=-h^2ω_K^2/f_π^2(e^-μ_Kr/4π r-e^-Λ' r/4π r-(Λ'^2-μ_K^2)e^-Λ' r/8πΛ')   ,

where we define Λ' as

    Λ'^2=Λ^2-q_0^2=Λ^2-ω_K^2   .

Using the above DD_s0^* potential and the DD potential provided
by the OBE model, we can check whether the three-body DDD_s0^*
system binds.

The binding energies we obtain with different cutoffs are tabulated
in Table <ref>.


With the effective cutoff Λ' ranging from 0.8-1.6 GeV,
the results of Table.<ref> indicate that the DDD_s0^*
bound state is located about (65-90)  MeV below
the DDDK threshold.

This is to be compared with 100  MeV for the full four-body calculation,
see Table.<ref> for details.

That is, as happened with the DD_s0^* / DDK system, the approximation
that the D_s0^* is a compact state results in underbinding
for the DDD_s0^* / DDDK system, but not much.



§ SUMMARY



In this manuscript we argued that the DK interaction is attractive
enough as to generate DK, DDK and DDDK bound states.

For this we began by assuming that the D_s0^*(2317) is a DK molecule,
which determines in turn the DK interaction.

Then, by means of the Gaussian Expansion Method <cit.> (a method for
few-body calculations),
we have addressed the question of whether one can build up multi-component
molecular states, similar to the formation of atomic nuclei
from clusters of nucleons bound by
the nucleon-nucleon interaction.

The answer is yes. We find a bound DDK trimer and a DDDK tetramer.

The prediction of this trimer confirms the previous calculations
of Refs. <cit.>,
while the prediction of the tetramer is novel
to the present work.

We have checked the robustness of these predictions
against a series of uncertainties.

While the DK interaction is well constrained by the existence of
the D_s0^*(2317) and chiral perturbation theory,
the DD interaction is considerably
less well-known.

Yet it also enters the calculations. We chose to describe the DD potential
in terms of the OBE model, in which the DD interaction turns out
to be mildly attractive and has a minor impact on the binding
energy of the trimer and tetramer states.

The DK potential, though well-known, is still subject to subleading
corrections, which we take into account by varying
the exact form of this potential.

As expected from the fact that we are dealing with subleading corrections,
the predictions are almost left unchanged by these variations.


In addition, we have studied a rather unlikely scenario that the D_s0^*(2317) is dominanty a genuine cs̅
state. Nonetheless, even in such a case, we still predict D D_s0^* and DDD_s0^* bound states
with the same quantum numbers as the DDK trimer and DDDK tetramer,
but this time located at approximately (50-62) and (60-90)  MeV
below the DDK and DDDK thresholds (instead of 70 and 100  MeV
when the D_s0^* is a molecular meson).

The binding mechanism is the long-range one-kaon-exchange potential
in the D D_s0^* system: owing to the mass difference between the
D and D_s0^* mesons, the kaon is exchanged near the mass shell,
leading to an enhancement in the range of
the potential <cit.>.

Although the existence of the DDK and DDDK bound states
seems to be quite robust, the question of where to find
them is much more challenging.

If we now focus on the DDK state, the experimental discovery of
the D_s0^*(2317) gives a clue.

As already argued in Ref. <cit.>,
but awaiting for a concrete study, the DDK state
can decay into D D_s^* or D^* D_s in P-wave.

Therefore one may look for inclusive combinations of three particles
D D_s π and search for structures in the corresponding invariant
mass distributions.

Given enough statistics, there should be a possibility to discover it
in the e^+ e^- collision data collected by Belle or BelleII or
in the pp collision data collected at the LHC. 

It is well known that heavy quark spin and flavor symmetries relate
the DK interaction to those of D^*K, BK̅ and B^*K̅.

This is consistent with the existence of the D_s1(2460).

The bottom counterparts of the D_s0^*(2317) and D_s1(2460)
have been predicted in a number of studies <cit.>
and confirmed by lattice QCD simulations <cit.>.

As a result, we naively expect the existence of the heavy quark
symmetry partners of the DDK and DDDK states.

At this moment, given the accessible center of mass energies at current
facilities, and the simplification that both the D and K 
are 0^- mesons that only decay weakly, we believe that they should
be of top priority both experimentally and theoretically.




§ ACKNOWLEDGEMENTS

This work is partly supported by the National Natural Science Foundation
of China under Grant No. 11735003, the Fundamental Research Funds
for the Central Universities, and the Thousand Talents Plan
for Young Professionals.



