# Molecular Interpretation of $X(3960)$ as $D_{s}^{+}D_{s}^{-}$ State

Halil Mutuk<EMAIL_ADDRESS>Department of Physics, Faculty of Sciences,
Ondokuz Mayis University, 55139 Samsun, Türkiye  Institute for Advanced
Simulation, Institut für Kernphysik and Jülich Center for Hadron Physics,
Forschungszentrum Jülich, D-52425 Jülich, Germany

###### Abstract

We study $D_{s}^{+}D_{s}^{-}$ and $D\bar{D}$ states assuming that they are
hadronic molecules with $J^{PC}=0^{++}$ quantum number. We use two-point QCD
sum rule formalism and extract the mass and decay constant values of these
states. We take into account contributions of various quark, gluon, and mixed
vacuum condensates up to dimension eight. The extracted mass and decay
constant values of $D\bar{D}$ and $D_{s}^{+}D_{s}^{-}$ states read as
$M_{D\bar{D}}=3795^{+85}_{-82}~{}\mathrm{MeV}$,
$f_{D\bar{D}}=1.70^{+0.33}_{-0.29}\times 10^{-2}~{}\mathrm{GeV}^{5}$, and
$M_{D_{s}^{+}D_{s}^{-}}=3983^{+93}_{-88}~{}\mathrm{MeV}$,
$f_{D_{s}^{+}D_{s}^{-}}=2.52^{+0.64}_{-0.54}\times
10^{-2}~{}\mathrm{GeV}^{5}$, respectively. The predicted mass of
$D_{s}^{+}D_{s}^{-}$ state is in good agreement with the recent LHCb
observation and supports quantum number and molecular picture assignments. A
possible observation of $D\bar{D}$ state would help for establishing the
lowest four-quark state in charmonium sector.

##  I Prologue

Multiquark states have been proposed when the quark model was introduced by
Gell-Mann [1] and Zweig [2] more than fifty years ago. After the proposition
of quark model, experimentally observed states fall in the scope of baryons
$(qqq)$ and mesons $(q\bar{q})$. Theoretical investigations which are
supported by the experimental observations of these traditional hadrons were
the success of the quark model. This situation turned into a different scene
when the observation of $\chi_{c1}(3872)$ (formerly known as $X(3872)$) state
announced by the Belle Collaboration in 2003 [3]. Apart from its mass which is
lower than the quark model prediction, this state decays into $J/\psi\rho$ and
$J/\psi\omega$ which are forbidden and OZI-suppressed for a bare $c\bar{c}$
configuration. Many more states are observed such as $Z_{c}(3900)$ [4],
$Z_{c}(4020)$ [5, 6], $P_{c}$ [7, 8] and $Z_{cs}(3985)$ [9] whose physical
properties cannot be obtained within the quark model framework.

Experimentally known charmonium spectrum lies in the range of 3.0-4.5 GeV. The
spectrum around 3.9 GeV would correspond to $\chi_{cJ}(2P)$ multiplet since
the spectrum of $\chi_{cJ}(1P)$ is in the range of 3.4-3.5 GeV and is quite
interesting because most of the experimentally observed charmoniumlike states
are near the $\bar{D}^{*0}D^{0}/\bar{D}^{0}D^{*0}$ threshold for
$\chi_{c1}(3872)$, $\bar{D}^{*}D/\bar{D}D^{*}$ threshold for $Z_{c}(3900)$,
$\bar{D}^{*}D^{*}$ threshold for $Z_{c}(4020)$ state, and
$\bar{D}_{s}D^{*}/\bar{D}_{s}^{*}D$ threshold for $Z_{cs}(3985)$ state.

In this respect, experimental observations put some puzzles through our
understanding of exotic states. In 2004 summer, the Belle Collaboration
observed a state with a mass around 3940 MeV in the invariant $\omega J/\psi$
mass distribution of the $B\to K\omega J/\psi$ [10]. The BaBar Collaboration
confirmed this state in the $B^{0,+}\to J/\psi\omega K^{0,+}$ but with a
different mass around 3915 MeV [11]. Determination of the quantum numbers for
this state came a couple of years later after the observation as
$J^{PC}=0^{++}$ by BaBar Collaboration [12]. In 2020, the LHCb Collaboration
found a structure named $\chi_{c0}(3915)$, in the $D^{+}D^{-}$ mass
distribution of $B^{+}\to D^{+}D^{-}K^{+}$ decay with the mass and width of
[13, 14]

$J^{PC}=0^{++};~{}M=3924\pm 2~{}\text{MeV};~{}\Gamma=17\pm 5~{}\text{MeV}.$
(1)

Particle Data Group refer these states as $\chi_{c0}(3915)$ [15]. This state
is attributed as conventional $\chi_{c0}(2P)$ charmonium state in Refs. [12,
16, 17]. Ref. [18] argued $\chi_{c0}(3915)$ as an $S$ wave $D_{s}\bar{D}_{s}$
molecular state with a binding energy of 18 MeV. It was mentioned in Ref. [19]
that experimental significance of the OZI forbidden $X(3915)\to\omega J/\psi$
decay could put a question mark to the conventional charmonia description of
this state.

Very recently, the LHCb Collaboration [20] announced the observation of a new
structure named as $X(3960)$ in the $D_{s}^{+}D_{s}^{-}$ mass distribution of
the $B^{+}\to D_{s}^{+}D_{s}^{-}K^{+}$ decay with the mass and width of

$J^{PC}=0^{++};~{}M=3955\pm 6\pm 11~{}\text{MeV};~{}\Gamma=48\pm 17\pm
11~{}\text{MeV}.$ (2)

The assigned quantum number makes this structure an $S$ wave state. The
resonant peak is just above the $D_{s}^{+}D_{s}^{-}$ threshold which is 3937
MeV and is a good candidate for a $D_{s}^{+}D_{s}^{-}$ molecule. The threshold
of $D\bar{D}$ is 3739 MeV. According to quark model calculations [21, 22],
$\chi_{c0}(2P)$ states have masses around 3.9 GeV and $\chi_{c0}(3P)$ states
have masses around 4.2 GeV. The quantum numbers, masses and widths of
$\chi_{c0}(3915)$ and $X(3960)$ suggest these particles as the same particle.
However, the branching ratio of these states

$\frac{\Gamma(X\to D^{+}D^{-})}{\Gamma(X\to D_{s}^{+}D_{s}^{-})}=0.29\pm
0.09\pm 0.10\pm 0.08$ (3)

says that $\Gamma(X\to D^{+}D^{-})<\Gamma(X\to D_{s}^{+}D_{s}^{-})$. This
implies the exotic nature of $X(3960)$ state because conventional charmonium
states predominantly decay into $D\bar{D}$ and $D^{*}\bar{D}^{*}$ states.
Furthermore it is harder to excite an $s\bar{s}$ pair from vacuum compared
with $u\bar{u}$ or $d\bar{d}$ pairs. Therefore, it is natural to expect that
$\Gamma(X\to D^{+}D^{-})$ should be bigger than $\Gamma(X\to
D_{s}^{+}D_{s}^{-})$.

The $D\bar{D}$ and $D_{s}^{+}D_{s}^{-}$ molecular states have been studied in
various approaches before [23, 24, 25, 26, 27, 28] and after [29, 30, 31, 32,
33, 34] the observation of $X(3960)$. In Refs. [23, 24, 25], a $D\bar{D}$
bound state was found whereas no bound state was found close to the
$D_{s}^{+}D_{s}^{-}$ threshold. If the $\chi_{c0}(3915)$ state couple both the
$D^{+}D^{-}$ and $D_{s}^{+}D_{s}^{-}$, it can produce an enhancement close to
the $D_{s}^{+}D_{s}^{-}$ threshold and would appear as a signal in the
$D_{s}^{+}D_{s}^{-}$ mass distribution. In that case, the experimental
observation of $X(3960)$ would be explained without introducing an extra
resonance. Refs. [26, 29] support this conclusion. In Ref. [30], using
nonrelativistic effective field theory it was found that a bound state or a
virtual state below the $D_{s}^{+}D_{s}^{-}$ mass threshold is needed to
describe the $D_{s}^{+}D_{s}^{-}$ mass distribution. Ref. [31] studied
$X(3960)$ as $D_{s}^{+}D_{s}^{-}$ molecule with $J^{PC}=0^{++}$ quantum number
in QCD sum rule framework. Production rates of $D\bar{D}$ and
$D_{s}^{+}D_{s}^{-}$ states were studied via effective Lagrangian approach
[32]. Using improved chromomagnetic interaction model, $S-$wave
$Qq\bar{Q}\bar{q}$ tetraquark states were studied in Ref. [33]. They observed
that $X(3960)$ state can be well explained in their model. In Ref. [34] it was
found that $\chi_{c0}(2P)$ state below the $D_{s}\bar{D}_{s}$ threshold can
depict the structure of $X(3960)$. In one-boson-exchange model
$D_{s}\bar{D}_{s}/D^{*}\bar{D}^{*}/D_{s}^{*}\bar{D}_{s}^{*}$ states were
studied considering the $S-D$ wave mixing effects [35].

In this present study, we calculate the mass and decay constant of the
$D^{+}D^{-}$ and $D_{s}^{+}D_{s}^{-}$ states in the context of the two-point
QCD sum rule (QCDSR) [36, 37]. We choose scalar interpolating current in
correlation function to account the quantum numbers of the $\chi_{c0}(3915)$
and $X(3960)$ state. We evaluate the two-point correlation function in terms
of various quark, gluon and mixed vacuum condensates.

The paper is organized in the following way: in Section II we provide the
method of QCDSR calculations for the considered states to extract mass and
coupling constant values. Section III presents the details of numerical
analyses and results for spectroscopic parameters. Section IV is reserved for
summary and conclusion.

##  II Two-Point Correlator

In order to extract mass and decay constant from QCDSR technique, an
interpolating current representing the physical state in the problem should be
constructed. To this end, we use the following interpolating currents
interpreting $D\bar{D}$ and $D_{s}^{+}D_{s}^{-}$ states as scalar molecules
with quantum number $J^{P}=0^{++}$

$\displaystyle j_{1}(x)$ $\displaystyle=$
$\displaystyle\bar{q}^{a}(x)\gamma_{\mu}c^{a}(x)\bar{c}^{b}\gamma^{\mu}q(x)^{b},$
(4) $\displaystyle j_{2}(x)$ $\displaystyle=$
$\displaystyle\bar{s}^{a}(x)\gamma_{\mu}c^{a}(x)\bar{c}^{b}\gamma^{\mu}s(x)^{b},$
(5)

where $q=u,d$, The sum rules for masses and decay constants of the $D\bar{D}$
and $D_{s}^{+}D_{s}^{-}$ states related to these currents are constructed from
the two-point correlation function

$\Pi(q)=i\int d^{4}xe^{iqx}\langle 0|T[j(x)j^{\dagger}(0)]|0\rangle,$ (6)

where $j(x)$ is the interpolating current, $T$ is the time ordered product
between the operators and $\langle\cdots\rangle$ is the QCD vacuum expectation
value.

The basic idea of the QCDSR framework is that there is an interval in momentum
$q$ where the correlation function may be equivalently described at both
quark-gluon (QCD or OPE) and hadron (phenomenological) levels. A typical QCDSR
calculation has two steps. In the first step, correlation function in the
phenomenological side is represented by the hadronic degrees of freedom.
Inserting complete set of intermediate states for the hadron of interest into
the correlation function Eq. (6) and isolating ground state contribution of
higher states and continuum, one can obtain

$\Pi^{\text{Phen}}(q)=\frac{\langle 0|j|H(p)\rangle\langle
H(p)|j^{\dagger}|0\rangle}{m_{H}^{2}-q^{2}}+\cdots.$ (7)

Let $m_{H}$ and $f_{H}$ be the mass and decay constant of the hadron,
respectively. One can write the matrix element as

$\langle 0|j|H\rangle=m_{H}f_{H}.$ (8)

With this definition $\Pi^{\text{Phen}}(q)$ as can be written as

$\Pi^{\text{Phen}}(q)=\frac{m_{H}^{2}f_{H}^{2}}{m_{H}^{2}-q^{2}}+\cdots,$ (9)

where $\cdots$ represents contributions of higher states. As a result of
writing the spectral function in terms of intermediate states for the hadron
of interest, all these states that couple to the interpolating current $j(x)$
contribute to the correlation function. In the case of four-quark mesons, it
should be mentioned that two-meson scattering states can contribute to the
$\Pi^{\text{Phen}}(q)$ in Eq. (7). Since these contributions are small, they
can be neglected [38, 39, 40, 41, 42, 43].

In the second step, correlation function in Eq. (6) is calculated via the
Operator Product Expansion (OPE) formulated in Ref. [44]. The OPE
representation of the correlation function $\Pi(q)$ can be obtained by
inserting interpolating current into the correlation function and contracting
heavy and light quark fields. This contraction is called Wick contraction. The
result is composed of heavy and light quark propagators:

$\displaystyle\Pi^{\text{OPE}}(q)$ $\displaystyle=$ $\displaystyle\int
d^{4}xe^{iqx}\text{Tr}[S^{a^{\prime}a}q_{(}-x)\gamma_{\mu}S^{aa^{\prime}}c(x)\gamma_{\nu}]$
(10) $\displaystyle\times$
$\displaystyle\text{Tr}[S^{b^{\prime}b}_{c}(-x)\gamma_{\mu}S^{bb^{\prime}}_{q}(x)\gamma_{\nu}].$

Here $S^{ab}_{q}$ and $S^{ab}_{c}$ are the light and heavy quark propagators,
respectively. The light quark propagator reads as

$\displaystyle S^{ab}_{q}(x)$ $\displaystyle=$ $\displaystyle
i\delta_{ab}\frac{x\\!\\!\\!/x}{2\pi^{2}x^{4}}-\delta_{ab}\frac{\langle\bar{q}q\rangle}{12}+i\delta_{ab}\frac{x\\!\\!\\!/m_{q}\langle\bar{q}q\rangle}{48}$
(11) $\displaystyle-$
$\displaystyle\delta_{ab}\frac{x^{2}}{192}\langle\bar{q}g_{s}\sigma
Gq\rangle+i\delta_{ab}\frac{x^{2}x\\!\\!\\!/m_{q}}{1152}\langle\bar{q}g_{s}\sigma
Gq\rangle$ $\displaystyle-$ $\displaystyle
ig_{s}\frac{G^{\alpha\beta}_{ab}}{32\pi^{2}x^{2}}[x\\!\\!\\!/\sigma_{\alpha\beta}+\sigma_{\alpha\beta}x\\!\\!\\!/]-i\delta_{ab}\frac{x^{2}x\\!\\!\\!/g_{s}^{2}\langle\bar{q}q\rangle^{2}}{7776}$
$\displaystyle-$
$\displaystyle\delta_{ab}x^{4}\frac{\langle\bar{q}q\rangle\langle
g_{s}^{2}G^{2}\rangle}{27648}+\cdots.$

In the light quark propagator we set $m_{q}\to 0$ for $q=u,d$. For the heavy
quark propagator $S^{ab}_{c}$, we employ the following expression

$\displaystyle S^{ab}_{c}(x)$ $\displaystyle=$ $\displaystyle
i\int\frac{d^{4}k}{(2\pi)^{4}}e^{-ikx}[\frac{\delta_{ab}(k\\!\\!\\!/+m_{c})}{k^{2}-m_{c}^{2}}$
(12) $\displaystyle-$
$\displaystyle\frac{g_{s}G^{\alpha\beta}_{ab}}{4}\frac{\sigma_{\alpha\beta}(k\\!\\!\\!/+m_{c})+(k\\!\\!\\!/+m_{c})\sigma_{\alpha\beta}}{(k^{2}-m_{c}^{2})^{2}}$
$\displaystyle+$
$\displaystyle\frac{g_{s}^{2}G^{2}}{12}\delta_{ab}m_{c}\frac{k^{2}+m_{c}k\\!\\!\\!/}{(k^{2}-m_{c}^{2})^{4}}+\cdots.]$

In Eqs. (11) and (12), we use the notations
$G^{\alpha\beta}_{ab}=G^{\alpha\beta}_{A}t^{A}_{\alpha\beta},~{}G^{2}=G^{A}_{\alpha\beta}G^{A}_{\alpha\beta}$,
where $a,b=1,2,3$ and $A,B,C=1,2,\cdots,8$ are color indices.
$t^{A}=\lambda^{A}/2$, and $\lambda^{A}$ are the Gell-Mann matrices. In the
nonperturbative terms the gluon field strength tensor
$G^{A}_{\alpha\beta}=G^{A}_{\alpha\beta}(0)$ is fixed at $x=0$.

The sum rules for mass $m_{H}$ and decay constant $f_{H}$ can be obtained via
equating the invariant amplitudes $\Pi^{\text{Phen}}(q)$ and
$\Pi^{\text{OPE}}(q)$; $\Pi^{\text{Phen}}(q)=\Pi^{\text{OPE}}(q)$. Such a
direct matching is weakened due to the OPE is valid at large $q^{2}$ and
phenomenological side is valid at small $q^{2}$. Furthermore, contributions of
higher states and continuum still exist. Borel transformation should be
applied in both sides of this equality in order to suppress contributions of
higher states and continuum. After this operation, there is still a remnant in
the sum rules. Using quark-hadron duality assumption, one subtracts higher
states and continuum terms from the phenomenological side of the matching (the
$\cdots$ terms in Eq. 9). After these standard manipulations of the QCDSR
method, the QPE side of the obtained sum rule takes the following form

$\tilde{\Pi}^{\text{QCD}}=\int_{\mathcal{M}^{2}}^{s_{0}}\rho(s)e^{-s/M^{2}}ds,$
(13)

where $\mathcal{M}=(2m_{c}+2m_{s})$ for $D_{s}^{+}D_{s}^{-}$,
$\mathcal{M}=(2m_{c})$ for $D\bar{D}$ (we neglect the mass of $u$ and $d$
quarks), $s_{0}$ is the continuum threshold which is the energy characterizing
beginning of the continuum, $M^{2}$ is Borel parameter and $\rho(s)$ is the
spectral density obtained from the imaginary part of the correlation function,
$\rho(s)=\text{Im}[\Pi(s)]/\pi$. We have taken into account vacuum condensates
up to dimension eight. The expressions for $\rho(s)$ are rather lengthy, so we
refrain to provide it here. Interested readers can request the spectral
densities from the author.

The phenomenological side of the sum rule has the following form

$\tilde{\Pi}^{\text{Phen}}=f_{H}^{2}m_{H}^{2}e^{-\frac{m_{H}^{2}}{M^{2}}}.$
(14)

Equating Eqs. (13) and (14), one can obtain sum rules for mass and decay
constant. The mass formula reads as

$m_{H}^{2}(s_{0},M^{2})=\frac{\tilde{\Pi}^{\prime\text{QCD}}}{\tilde{\Pi}^{\text{QCD}}}$
(15)

where
$\tilde{\Pi}^{\prime\text{QCD}}=\frac{d}{d(-\frac{1}{M^{2}})}\tilde{\Pi}^{\text{QCD}}$,
and decay constant formula reads as

$f_{H}^{2}(s_{0},M^{2})=\frac{e^{m_{H}^{2}/M^{2}}}{m_{H}^{2}}\tilde{\Pi}^{\text{QCD}}.$
(16)

##  III Spectroscopic Parameters of $D\bar{D}$ and $D_{s}^{+}D_{s}^{-}$ States

There are two parameter groups for a QCDSR calculation; QCD input parameters
and resulting sum rule parameters of the model. The first group consists of
quark masses and QCD condensates which are given as

$\displaystyle\langle\overline{q}q\rangle=-(0.24\pm
0.01)^{3}~{}\mathrm{GeV}^{3},\ \langle\overline{s}s\rangle=(0.8\pm
0.1)\langle\overline{q}q\rangle,$ $\displaystyle\langle\overline{q}g_{s}\sigma
Gq\rangle=m_{0}^{2}\langle\overline{q}q\rangle,\
\langle\overline{s}g_{s}\sigma
Gs\rangle=m_{0}^{2}\langle\overline{s}s\rangle,$ $\displaystyle
m_{0}^{2}=(0.8\pm 0.2)~{}\mathrm{GeV}^{2}$
$\displaystyle\langle\frac{\alpha_{s}G^{2}}{\pi}\rangle=(0.012\pm
0.004)~{}\mathrm{GeV}^{4},$ $\displaystyle
m_{s}=93_{-5}^{+11}~{}\mathrm{MeV},\ m_{c}=1.27\pm 0.2~{}\mathrm{GeV}.$ (17)

The second group contains the continuum threshold parameter $s_{0}$ which
arises after continuum subtraction via using quark-hadron duality and the
$M^{2}$ Borel parameter which arises after Borel transformation. These are
auxiliary parameters and any quantity which is extracted from QCDSR
calculations should be independent or have a mild dependence with respect to
variation of these parameters. Notwithstanding, these auxiliary parameters
should meet some criteria of the QCDSR method. Since there is no single value
for $s_{0}$ and $M^{2}$, a suitable working regions of these parameters must
be found in order to get reliable predictions. Finding proper choice of
continuum threshold parameter $s_{0}$ and Borel parameter $M^{2}$ is an
important task in the sum rule analysis.

The region for Borel parameter $M^{2}$ is determined by OPE convergence and
suppression of the contributions arising from the higher and continuum states.
Finding upper and lower limits of Borel parameter requires additional
elaboration of the sum rules. Pole contribution (PC),

$\text{PC}=\frac{\Pi^{\text{OPE}}(s_{0},M^{2})}{\Pi^{\text{OPE}}(s_{0}=\infty,M^{2})},$
(18)

and convergence of the operator product expansion

$R(M^{2})=\frac{\Pi^{\text{OPE,DimN}}(s_{0},M^{2})}{\Pi^{\text{OPE}}(s_{0}=\infty,M^{2})}.$
(19)

will help to determine upper and lower limits of the Borel parameter. In Eq.
19, $\Pi^{\text{OPE,DimN}}(s_{0},M^{2})$ is the expression used to estimate
convergence of OPE. In QCD sum rule approach, $\mathcal{O}(\alpha_{s})$
corrections to leading and next-to-leading order terms in
$\Pi^{\text{OPE}}(s_{0},M^{2})$ maybe comparable above Dim 8 dimensional
contributions. Due to the interpolating current we used, contribution of Dim8
disappears due to the gamma matrix algebra. Therefore we use
$\Pi^{\text{Dim7}}(s_{0},M^{2})$ to look for the convergence of OPE. In this
way, we can avoid uncertainties connected to $\mathcal{O}(\alpha_{s})$
corrections.

For the upper limit, PC should reside fairly in the correlation function while
for the lower limit PC must be dominant. This is due to the fact that after
Borel transformation, an exponential term emerges in both side of the sum
rule, as can be seen from Eqs. 13 and 14. The integrals receive main
contribution at $s\simeq M^{2}$. For lower $M^{2}$ values, the integral (Eq.
13) receives contributions from lower $s$ values and for higher $M^{2}$ values
these contributions decrease. $R(M^{2})$ can be used to determine the lower
limit of Borel parameter. Between upper and lower limits of $M^{2}$ a balance
can be maintained where dominance of single resonance contribution and
convergence of the series take place.

The continuum threshold $s_{0}$ eliminates contributions from the continuum
states and higher resonances, and isolates the ground state in the correlation
function. $s_{0}$ is the place where the continuum states and higher
resonances begin to contribute to the correlation function. Therefore $s_{0}$
specifies the threshold in the integral which needs to be carried out. Since
in the sum rule calculations ground states are being handled, continuum
threshold $s_{0}$ starts with the first excited state. The energy that needs
to excite the ground state is $\delta=\sqrt{s_{0}}-m$ where $m$ is the mass of
the ground state. In conventional mesons and baryons, parameters of first
excited states are known experimentally or have a theoretical consensus. Based
on this knowledge, it is quite possible to estimate $s_{0}$ value. For
experimentally known mesons and baryons $\delta$ varies in the range of
$0.3~{}\text{GeV}\leq\delta\leq 0.8~{}\text{GeV}$. In the case of multiquark
hadrons (exotic states) there can be lack of relevant information for excited
states and ground states. The determination of $s_{0}$ can be obtained via
aforementioned limits on PC and the convergence of OPE. The region for
continuum threshold $s_{0}$ can be chosen as possible as small to provide a
reliable Borel region and meanwhile to obtain a maximum for PC. With all these
considerations it is possible to examine obtained sum rules numerically. Our
analysis gives

$\displaystyle M^{2}\in[2.0,2.4]~{}\mathrm{GeV}^{2},\
s_{0}\in[17,19]~{}\mathrm{GeV}^{2},$ (20) $\displaystyle
M^{2}\in[2.3,2.7]~{}\mathrm{GeV}^{2},\ s_{0}\in[19,21]~{}\mathrm{GeV}^{2},$
(21)

for $D\bar{D}$ and $D_{s}^{+}D_{s}^{-}$, respectively. The working regions of
Eqs. (20) and (21) meet all the restrictions mentioned before. A further
discussion should be done for the continuum threshold $s_{0}$ value.
Considering $D^{*\pm}(2010)$ state for the first excited state of $D$, the
mass [45] square of two-$D^{*\pm}(2010)$ meson gives
$(2m_{D^{*\pm}})^{2}\simeq 16.16~{}\mathrm{GeV}^{2}$. On the otherside,
considering $D^{*0}(2007)$ state for the first excited state of $D$, the mass
[45] square of two-$D^{*0}(2007)$ meson gives $(2m_{D^{*0}})^{2}\simeq
16.10~{}\mathrm{GeV}^{2}$. By same token, considering $D_{s}^{*\pm}$ state as
the first excited state of $D_{s}^{\pm}$, the mass [45] square of
two-$D_{s}^{*\pm}$ meson gives $(2m_{D_{s}^{*\pm}})^{2}\simeq
17.84~{}\mathrm{GeV}^{2}$. The obtained continuum threshold $s_{0}$ values are
reasonable to excite $D\bar{D}$ and $D_{s}^{+}D_{s}^{-}$ to their first
excited states.

Working regions of Borel parameter $M^{2}$ and continuum threshold $s_{0}$
give the pole contributions as

$\displaystyle D\bar{D}$ $\displaystyle:$ $\displaystyle 0.32\leq\text{PC}\leq
0.62,$ (22) $\displaystyle D_{s}^{+}D_{s}^{-}$ $\displaystyle:$ $\displaystyle
0.29\leq\text{PC}\leq 0.73,$ (23)

for upper and lower limits of Borel parameter $M^{2}$, respectively. In a
standard QCDSR analysis, the minimum value for PC of four-quark systems is
approximately 0.15-0.2 [46]. At the upper limit of Borel parameter $M^{2}$, PC
should exceed this range. After providing PC in an appropriate range, one can
extract the mass and decay constant values for the corresponding state. Our
results for masses and decay constants of $D\bar{D}$ and $D_{s}^{+}D_{s}^{-}$
states read as

$\displaystyle M_{D\bar{D}}$ $\displaystyle=$ $\displaystyle
3795^{+85}_{-82}~{}\mathrm{MeV},$ $\displaystyle f_{D\bar{D}}$
$\displaystyle=$ $\displaystyle 1.70^{+0.33}_{-0.29}\times
10^{-2}~{}\mathrm{GeV}^{5},$ (24)

and

$\displaystyle M_{D_{s}^{+}D_{s}^{-}}$ $\displaystyle=$ $\displaystyle
3983^{+93}_{-88}~{}\mathrm{MeV},$ $\displaystyle f_{D_{s}^{+}D_{s}^{-}}$
$\displaystyle=$ $\displaystyle 2.52^{+0.64}_{-0.54}\times
10^{-2}~{}\mathrm{GeV}^{5}.$ (25)

In Figures 1, 2, 3 and 4 we plot the mass and decay constant of $D\bar{D}$ and
$D_{s}^{+}D_{s}^{-}$ states as functions of $M^{2}$ and $s_{0}$, respectively.

Figure 1: The mass of the $D\bar{D}$ state as a function of $M^{2}$ at fixed
$s_{0}$ (left panel), and as a function of $s_{0}$ at fixed $M^{2}$ (right
panel).

Figure 2: Same as in Fig. 1, but for the coupling $f$ of the $D\bar{D}$ state.

Figure 3: The mass of the $D_{s}^{+}D_{s}^{-}$ state as a function of $M^{2}$
at fixed $s_{0}$ (left panel), and as a function of $s_{0}$ at fixed $M^{2}$
(right panel).

Figure 4: Same as in Fig. 3, but for the coupling $f$ of the
$D_{s}^{+}D_{s}^{-}$ state.

As can be seen from Figs. 1 and 3, the mass values are very stable with
respect to changes in Borel parameter $M^{2}$ and continuum threshold $s_{0}$.
Figs. 2 and 4 have a residual dependence on the Borel parameter $M^{2}$ and
continuum threshold $s_{0}$. This dependence is an essential part of the
theoretical errors in QCD sum rule calculations. The reason for this situation
could be as a result of mass and decay constant extractions from the QCD sum
rules. The mass of the state is determined as
$m^{2}=\frac{\Pi^{\prime}(s_{0},M^{2})}{\Pi(s_{0},M^{2})}$. The spectral
densities in both numerator and denominator could smooth variations of $M^{2}$
and $s_{0}$. On the other side, the decay constant is determined as
$f^{2}=\frac{e^{m^{2}/M^{2}}}{m^{2}}\Pi(s_{0},M^{2})$ where the numerical
value can change with respect to changes in $M^{2}$ and $s_{0}$ since such a
ratio (as in the case of mass) of correlation functions is absent. Theoretical
uncertainties of the extracted masses are 2.1 % for $D\bar{D}$ state and 2.1 %
for $D_{s}^{+}D_{s}^{-}$ state. In the decay constant values, theoretical
uncertainties are 18 % for $D\bar{D}$ state and 23 % for $D_{s}^{+}D_{s}^{-}$
state. These theoretical uncertainties are well below the accepted limits in
QCDSR calculations.

The extracted mass $M_{D_{s}^{+}D_{s}^{-}}=3983^{+93}_{-88}~{}\mathrm{MeV}$
for $D_{s}^{+}D_{s}^{-}$ state with the $J^{PC}=0^{++}$ quantum number agree
well with the experimental mass $M=3955\pm 6\pm 11~{}\text{MeV}$ of $X(3960)$.
Our prediction supports $X(3960)$ to be a $D_{s}^{+}D_{s}^{-}$ molecular state
with $J^{PC}=0^{++}$ quantum number. A recent QCDSR study obtained mass of
$D_{s}^{+}D_{s}^{-}$ state as $M=3.98\pm 0.10~{}\mathrm{GeV}$ and decay
constant $f=2.36^{+0.45}_{-0.45}\times 10^{-2}~{}\mathrm{GeV}^{5}$, with
$J^{PC}=0^{++}$ in the molecular picture [31]. Our results agree well with
their predictions.

The extracted mass for $D\bar{D}$ state with the $J^{PC}=0^{++}$ quantum
number is $M_{D\bar{D}}=3795^{+85}_{-82}~{}\mathrm{MeV}$. Ref. [31] obtained
mass and decay constant of $D\bar{D}$ state as $M=3.74\pm 0.09~{}\mathrm{GeV}$
and $f=1.61^{+0.23}_{-0.23}\times 10^{-2}~{}\mathrm{GeV}^{5}$. Our results
agree well with their predictions. The existence of such a bound state near
the $D\bar{D}$ threshold was demonstrated by analysis of $\gamma\gamma\to
D\bar{D}$ [47, 48]. The extracted mass of $D\bar{D}$ state lies nearby by the
$\psi(3770)$ with mass $M=3773.7\pm 0.4~{}\mathrm{MeV}$. The quantum number of
$\psi(3770)$ is $J^{PC}=1^{++}$. Conventional $\chi_{c0}(1P)$ state has a mass
of $M=3414.71\pm 0.30~{}\mathrm{MeV}$ with $J^{PC}=0^{++}$ quantum number. The
locations of these states can be seen in Fig. 5. The $D\bar{D}$ state has no
overlaps with these states. Therefore possible observation of $D\bar{D}$ state
is important in terms of establishing the lowest four-quark states in
charmonium sector.

Figure 5:  Predicted mass values of $D\bar{D}$ and $D_{s}^{+}D_{s}^{-}$
states, and the lowest $\chi_{c0}(1P)$ state together with $\psi(3770)$. The
lowest state in $J^{PC}=0^{++}$ is $\chi_{c0}(1P)$ and the state in
$J^{PC}=1^{++}$ is $\psi(3770)$. 3795 and 3983 refer to mass values of
$D\bar{D}$ and $D_{s}^{+}D_{s}^{-}$ states in $J^{PC}=0^{++}$, respectively.

##  IV Epilogue

In this present study, we considered $D\bar{D}$ and $D_{s}^{+}D_{s}^{-}$
states with $J^{PC}=0^{++}$ quantum number in molecular picture. We extracted
mass and decay constant (residue) values by using two-point QCDSR formalism.
We carried out the OPE up to the vacuum condensates of dimension eight.
Establishing pole contribution within the limits of QCDSR, we determined
working regions of continuum threshold $s_{0}$ and Borel parameter $M^{2}$.
The extracted mass and decay constant values of $D\bar{D}$ and
$D_{s}^{+}D_{s}^{-}$ states read as
$M_{D\bar{D}}=3795^{+85}_{-82}~{}\mathrm{MeV}$,
$f_{D\bar{D}}=1.70^{+0.33}_{-0.29}\times 10^{-2}~{}\mathrm{GeV}^{5}$, and
$M_{D_{s}^{+}D_{s}^{-}}=3983^{+93}_{-88}~{}\mathrm{MeV}$,
$f_{D_{s}^{+}D_{s}^{-}}=2.52^{+0.64}_{-0.54}\times
10^{-2}~{}\mathrm{GeV}^{5}$, respectively. The predicted mass for
$D_{s}^{+}D_{s}^{-}$ state agrees well with the experimental value $M=3955\pm
6\pm 11~{}\text{MeV}$ of the very recent LHCb observation of $X(3960)$. Our
result for $D_{s}^{+}D_{s}^{-}$ state supports the quantum number assignment
and molecular picture of $X(3960)$. On the other side, the possible
observation of $D\bar{D}$ state may pave the fay for the lowest four-quark
state in charmonium sector. We hope our outcomes may be helpful for future
theoretical studies and experimental endeavours.

###### Acknowledgements.

The author of this work is supported by The Scientific and Technological
Research Council of Türkiye (TUBITAK) in the framework of BIDEB-2219
International Postdoctoral Research Fellowship Program.

## References

  * Gell-Mann [1964] M. Gell-Mann, Phys. Lett. 8, 214 (1964).
  * Zweig [1964] G. Zweig, (1964).
  * Choi  _et al._ [2003] S. K. Choi _et al._ (Belle), Phys. Rev. Lett. 91, 262001 (2003), arXiv:hep-ex/0309032 .
  * Bondar  _et al._ [2012] A. Bondar _et al._ (Belle), Phys. Rev. Lett. 108, 122001 (2012), arXiv:1110.2251 [hep-ex] .
  * Ablikim  _et al._ [2013] M. Ablikim _et al._ (BESIII), Phys. Rev. Lett. 111, 242001 (2013), arXiv:1309.1896 [hep-ex] .
  * Ablikim  _et al._ [2014] M. Ablikim _et al._ (BESIII), Phys. Rev. Lett. 112, 132001 (2014), arXiv:1308.2760 [hep-ex] .
  * Aaij  _et al._ [2015] R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 115, 072001 (2015), arXiv:1507.03414 [hep-ex] .
  * Aaij  _et al._ [2019] R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 122, 222001 (2019), arXiv:1904.03947 [hep-ex] .
  * Ablikim  _et al._ [2021] M. Ablikim _et al._ (BESIII), Phys. Rev. Lett. 126, 102001 (2021), arXiv:2011.07855 [hep-ex] .
  * Abe  _et al._ [2005] K. Abe _et al._ (Belle), Phys. Rev. Lett. 94, 182002 (2005), arXiv:hep-ex/0408126 .
  * Aubert  _et al._ [2008] B. Aubert _et al._ (BaBar), Phys. Rev. Lett. 101, 082001 (2008), arXiv:0711.2047 [hep-ex] .
  * Lees  _et al._ [2012] J. P. Lees _et al._ (BaBar), Phys. Rev. D 86, 072002 (2012), arXiv:1207.2651 [hep-ex] .
  * Aaij  _et al._ [2020a] R. Aaij _et al._ (LHCb), Phys. Rev. Lett. 125, 242001 (2020a), arXiv:2009.00025 [hep-ex] .
  * Aaij  _et al._ [2020b] R. Aaij _et al._ (LHCb), Phys. Rev. D 102, 112003 (2020b), arXiv:2009.00026 [hep-ex] .
  * Zyla  _et al._ [2020] P. A. Zyla _et al._ (Particle Data Group), PTEP 2020, 083C01 (2020).
  * Duan  _et al._ [2020] M.-X. Duan, S.-Q. Luo, X. Liu, and T. Matsuki, Phys. Rev. D 101, 054029 (2020), arXiv:2002.03311 [hep-ph] .
  * Duan  _et al._ [2021] M.-X. Duan, J.-Z. Wang, Y.-S. Li, and X. Liu, Phys. Rev. D 104, 034035 (2021), arXiv:2104.09132 [hep-ph] .
  * Li and Voloshin [2015] X. Li and M. B. Voloshin, Phys. Rev. D 91, 114014 (2015), arXiv:1503.04431 [hep-ph] .
  * González [2017] P. González, J. Phys. G 44, 075004 (2017), arXiv:1611.03710 [hep-ph] .
  * E. Spadaro Norella and C. Chen (2022) [LHCb] E. Spadaro Norella and C. Chen (LHCb), “Particle zoo 2.0: New tetra- and pentaquarks at lhcb,” url: https://indico.cern.ch/event/1176505 (2022).
  * Barnes  _et al._ [2005] T. Barnes, S. Godfrey, and E. S. Swanson, Phys. Rev. D 72, 054026 (2005), arXiv:hep-ph/0505002 .
  * Li and Chao [2009] B.-Q. Li and K.-T. Chao, Phys. Rev. D 79, 094004 (2009), arXiv:0903.5506 [hep-ph] .
  * Gamermann  _et al._ [2007] D. Gamermann, E. Oset, D. Strottman, and M. J. Vicente Vacas, Phys. Rev. D 76, 074016 (2007), arXiv:hep-ph/0612179 .
  * Nieves and Valderrama [2012] J. Nieves and M. P. Valderrama, Phys. Rev. D 86, 056004 (2012), arXiv:1204.2790 [hep-ph] .
  * Hidalgo-Duque  _et al._ [2013] C. Hidalgo-Duque, J. Nieves, and M. P. Valderrama, Phys. Rev. D 87, 076006 (2013), arXiv:1210.5431 [hep-ph] .
  * Prelovsek  _et al._ [2021] S. Prelovsek, S. Collins, D. Mohler, M. Padmanath, and S. Piemonte, JHEP 06, 035 (2021), arXiv:2011.02542 [hep-lat] .
  * Meng  _et al._ [2021] L. Meng, B. Wang, and S.-L. Zhu, Sci. Bull. 66, 1288 (2021), arXiv:2012.09813 [hep-ph] .
  * Dong  _et al._ [2021] X.-K. Dong, F.-K. Guo, and B.-S. Zou, Progr. Phys. 41, 65 (2021), arXiv:2101.01021 [hep-ph] .
  * Bayar  _et al._ [2022] M. Bayar, A. Feijoo, and E. Oset, (2022), arXiv:2207.08490 [hep-ph] .
  * Ji  _et al._ [2022] T. Ji, X.-K. Dong, M. Albaladejo, M.-L. Du, F.-K. Guo, and J. Nieves, Phys. Rev. D 106, 094002 (2022), arXiv:2207.08563 [hep-ph] .
  * Xin  _et al._ [2022] Q. Xin, Z.-G. Wang, and X.-S. Yang, (2022), arXiv:2207.09910 [hep-ph] .
  * Xie  _et al._ [2022] J.-M. Xie, M.-Z. Liu, and L.-S. Geng, (2022), arXiv:2207.12178 [hep-ph] .
  * Guo  _et al._ [2022a] T. Guo, J. Li, J. Zhao, and L. He, (2022a), arXiv:2211.10834 [hep-ph] .
  * Guo  _et al._ [2022b] D. Guo, J.-Z. Wang, D.-Y. Chen, and X. Liu, Phys. Rev. D 106, 094037 (2022b), arXiv:2210.16720 [hep-ph] .
  * Chen and Huang [2022] R. Chen and Q. Huang, (2022), arXiv:2209.05180 [hep-ph] .
  * Shifman  _et al._ [1979a] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B 147, 385 (1979a).
  * Shifman  _et al._ [1979b] M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B 147, 448 (1979b).
  * Lee  _et al._ [2005] S. H. Lee, H. Kim, and Y. Kwon, Phys. Lett. B 609, 252 (2005), arXiv:hep-ph/0411104 .
  * Wang [2015] Z.-G. Wang, Int. J. Mod. Phys. A 30, 1550168 (2015), arXiv:1502.01459 [hep-ph] .
  * Sundu  _et al._ [2019] H. Sundu, S. S. Agaev, and K. Azizi, Eur. Phys. J. C 79, 215 (2019), arXiv:1812.10094 [hep-ph] .
  * Agaev  _et al._ [2019] S. S. Agaev, K. Azizi, B. Barsbay, and H. Sundu, Nucl. Phys. B 939, 130 (2019), arXiv:1806.04447 [hep-ph] .
  * Wang [2020a] Z.-G. Wang, Phys. Rev. D 101, 074011 (2020a), arXiv:2001.04095 [hep-ph] .
  * Wang [2020b] Z.-G. Wang, (2020b), arXiv:2005.12735 [hep-ph] .
  * Wilson [1969] K. G. Wilson, Phys. Rev. 179, 1499 (1969).
  * Workman  _et al._ [2022] R. L. Workman _et al._ (Particle Data Group), PTEP 2022, 083C01 (2022).
  * Agaev  _et al._ [2020] S. Agaev, K. Azizi, and H. Sundu, Turk. J. Phys. 44, 95 (2020), arXiv:2004.12079 [hep-ph] .
  * Wang  _et al._ [2021] E. Wang, H.-S. Li, W.-H. Liang, and E. Oset, Phys. Rev. D 103, 054008 (2021), arXiv:2010.15431 [hep-ph] .
  * Deineka  _et al._ [2022] O. Deineka, I. Danilkin, and M. Vanderhaeghen, Phys. Lett. B 827, 136982 (2022), arXiv:2111.15033 [hep-ph] .