# Near-threshold structures in the $D_{s}^{+}D_{s}^{-}$ mass distribution of
the decay $B^{+}\rightarrow D_{s}^{+}D_{s}^{-}K^{+}$

S. S. Agaev  Institute for Physical Problems, Baku State University, Az–1148
Baku, Azerbaijan  K. Azizi  Department of Physics, University of Tehran, North
Karegar Avenue, Tehran 14395-547, Iran  Department of Physics, Doǧuş
University, Dudullu-Ümraniye, 34775 Istanbul, Turkey  H. Sundu  Department of
Physics Engineering, Istanbul Medeniyet University, 34700 Istanbul, Turkey
Department of Physics, Kocaeli University, 41380 Izmit, Turkey

###### Abstract

Two near-threshold peaking structures with spin-parities
$J^{\mathrm{PC}}=0^{++}$ were recently discovered by the LHCb collaboration in
the $D_{s}^{+}D_{s}^{-}$ invariant mass distribution of the decay
$B^{+}\rightarrow D_{s}^{+}D_{s}^{-}K^{+}$. The first of them is the resonance
$X(3960)$, whereas the second one $X_{0}(4140)$ is a structure with the mass
around $4140~{}\mathrm{MeV}$. To explore their natures and model them, we
study the hadronic molecule $\mathcal{M}=D_{s}^{+}D_{s}^{-}$ and calculate its
mass, current coupling, and width. The mass and current coupling of the
molecule are extracted from the QCD two-point sum rule analyses by taking into
account vacuum condensates up to dimension $10$. To evaluate its full width,
we consider the processes $\mathcal{M}\to D_{s}^{+}D_{s}^{-}$ and
$\mathcal{M}\to\eta_{c}\eta^{(\prime)}$. Partial widths of these decays are
determined by the strong couplings $g_{i},\ i=1,2,3$ at vertices
$\mathcal{M}D_{s}^{+}D_{s}^{-}$ and $\mathcal{M}\eta_{c}\eta^{(\prime)}$. They
are computed by means of the three-point sum rule method. Predictions for the
mass $m=(4117\pm 85)~{}\mathrm{MeV}$ and width $\Gamma_{\mathcal{M}}=(59\pm
12)~{}\mathrm{MeV}$ of the molecule $\mathcal{M}$ are compared with the
corresponding LHCb data, and also with our results for the diquark-antidiquark
state $X=[cs][\overline{c}\overline{s}]$. We argue that the structure
$X_{0}(4140)$ may be interpreted as the hadronic molecule
$D_{s}^{+}D_{s}^{-}$, whereas the resonance $X(3960)$ can be identified with
the tetraquark $X$.

##  I Introduction

Different $X$ resonances discovered and studied during last years by the LHCb
collaboration became unavoidable and important part of exotic hadron
spectroscopy. Thus, the resonances $X(4140)$, $X(4274)$, $X(4500)$ and
$X(4700)$ were seen in the process $B^{+}\rightarrow J/\psi\phi K^{+}$ as
peaks in the $J/\psi\phi$ invariant mass distribution Aaij:2016iza . The
structures $X(4140)$ and $X(4274)$ bear the quantum numbers
$J^{\mathrm{PC}}=1^{++}$, and in the tetraquark model are composed of the
quarks $c\overline{c}s\overline{s}$. The resonances $X(4500)$ and $X(4700)$
are scalar particles with spin-parities $J^{\mathrm{PC}}=0^{++}$ and the same
quark content. It is worth noting, that $X(4140)$ and $X(4274)$ were observed
previously in the decays $B^{\pm}\rightarrow J/\psi\phi K^{\pm}$ by different
collaborations Aaltonen:2009tz ; Chatrchyan:2013dma ; Abazov:2013xda , whereas
the scalar resonances $X(4500)$ and $X(4700)$ were discovered for the first
time by LHCb. The resonance $X(4630)$ fixed in the $J/\psi\phi$ invariant mass
distribution of the decay $B^{+}\rightarrow J/\psi\phi K^{+}$ is a vector
member of the $X$ tetraquarks’ family LHCb:2021uow .

Recently, LHCb reported about new hidden charm-strange structures in the
$D_{s}^{+}D_{s}^{-}$ invariant mass distribution of the decay
$B^{+}\rightarrow D_{s}^{+}D_{s}^{-}K^{+}$ LHCb:2022vsv . One of them
$X(3960)$ is presumably a tetraquark $c\overline{c}s\overline{s}$ with quantum
numbers $J^{\mathrm{PC}}=0^{++}$, and the following parameters

$\displaystyle m_{1\mathrm{exp}}$ $\displaystyle=$ $\displaystyle(3956\pm 5\pm
10)~{}\mathrm{MeV},$ $\displaystyle\Gamma_{1\mathrm{exp}}$ $\displaystyle=$
$\displaystyle(43\pm 13\pm 8)~{}\mathrm{MeV}.$ (1)

The LHCb also found evidence for a second structure $X_{0}(4140)$ with the
mass around $4140~{}\mathrm{MeV}$. The mass and full width of this state are

$\displaystyle m_{2\mathrm{exp}}$ $\displaystyle=$ $\displaystyle(4133\pm 6\pm
6)~{}\mathrm{MeV},$ $\displaystyle\Gamma_{2\mathrm{exp}}$ $\displaystyle=$
$\displaystyle(67\pm 17\pm 7)~{}\mathrm{MeV}.$ (2)

The $X_{0}(4140)$ may be interpreted either as a new resonance with
$J^{\mathrm{PC}}=0^{++}$ assignment or $J/\psi\phi\leftrightarrow
D_{s}^{+}D_{s}^{-}$ coupled-channel effect LHCb:2022vsv . In the present work,
we assume that the structure $X_{0}(4140)$ is a second resonance seen by LHCb
in the $D_{s}^{+}D_{s}^{-}$ mass distribution.

The $X$ resonances are interesting objects for theoretical investigations:
Features of exotic mesons $c\overline{c}s\overline{s}$ were studied in
numerous articles by employing different models and technical methods
Nieves:2012tt ; Wang:2013exa ; Lebed:2016yvr ; Chen:2017dpy ; Agaev:2017foq ;
Meng:2020cbk ; Sundu:2018toi ; Agaev:2022iha . Some of these states was
undergone to rather detailed exploration also in our publications. Thus, the
axial-vector resonances $X(4140)$ and $X(4274)$ were analyzed in Ref.
Agaev:2017foq , in which they were modeled as states composed of scalar and
axial-vector (anti)diquarks. In the case of $X(4140)$ the constituent diquark
(antidiquark) is the antitriplet (triplet) state of the color group
$SU_{c}(3)$, whereas to model $X(4274)$ we used (anti)diquarks from the sextet
representation of $SU_{c}(3)$. Predictions for masses and full widths of these
tetraquarks were compared with the LHCb data. It turned out, for masses of the
resonances $X(4140)$ and $X(4274)$ these models led to nice agreements with
the LHCb data. The model based on color triplet (anti)diquarks reproduced also
the full width of $X(4140)$, therefore it could be considered as a serious
candidate to resonance $X(4140)$. The width of the construction with color
sextet constituents is wider than that of $X(4274)$, which excludes it from a
list of possible pretenders.

The vector resonances $Y(4660)$ and $X(4630)$ were explored in our articles
[see, Refs. Sundu:2018toi ; Agaev:2022iha ] as tetraquarks
$[cs][\overline{c}\overline{s}]$ with spin-parities $J^{\mathrm{PC}}=1^{--}$
and $J^{\mathrm{PC}}=1^{-+}$, respectively. Predictions for the masses and
full widths of these states allowed us to interpret $Y(4660)$ and $X(4630)$ as
exotic mesons with diquark-antidiquark composition.

The discovery of structures $X(3960)$ and $X_{0}(4140)$ activated
investigations of hidden charm-strange resonances to account for features of
these new states Bayar:2022dqa ; Ji:2022uie ; Xin:2022bzt ; Xie:2022lyw ;
Chen:2022dad ; Guo:2022ggl . The $X(3960)$ was considered as a coupled-channel
effect Bayar:2022dqa , or as near the $D_{s}^{+}D_{s}^{-}$ threshold
enhancement by the conventional $P$-wave charmonium $\chi_{c0}(2P)$
Guo:2022ggl . The hadronic $D_{s}^{+}D_{s}^{-}$ molecule model was suggested
in Ref. Xin:2022bzt  to explain observed properties of the resonance
$X(3960)$.

In our paper Agaev:2022pis , we examined the tetraquark
$X=[cs][\overline{c}\overline{s}]$ with quantum numbers
$J^{\mathrm{PC}}=0^{++}$ and calculated its mass and full width. The
spectroscopic parameters of $X$, i.e., its mass and current coupling, were
found by means of the QCD two-point sum rule method. The full width of this
state was evaluated using decay channels $X\rightarrow D_{s}^{+}D_{s}^{-}$ and
$X\rightarrow\eta_{c}\eta^{(\prime)}$. Predictions for the mass $m=(3976\pm
85)~{}\mathrm{MeV}$ and width $\Gamma_{\mathrm{X}}=(42.2\pm
8.3)~{}\mathrm{MeV}$ obtained in Ref. Agaev:2022pis  allowed us to consider
the diquark-antidiquark state $X$ as an acceptable model for $X(3960)$.

In the present article, we continue our studies of the resonance $X(3960)$ and
include into analysis also the structure $X_{0}(4140)$. We investigate the
hadronic molecule $\mathcal{M}=D_{s}^{+}D_{s}^{-}$ by computing its mass,
current coupling, and full width. The mass and current coupling of this state
are calculated in the context of QCD two-point sum rule approach. The full
width of $\mathcal{M}$ is estimated by considering the decay channels
$\mathcal{M}\rightarrow D_{s}^{+}D_{s}^{-}$,
$\mathcal{M}\rightarrow\eta_{c}\eta$ and
$\mathcal{M}\rightarrow\eta_{c}\eta^{\prime}$. Partial widths of these
processes, apart from parameters of the initial and final particles, depend
also on strong couplings $g_{i},\ i=1,2,3$ at vertices
$\mathcal{M}D_{s}^{+}D_{s}^{-}$, $\mathcal{M}\eta_{c}\eta^{\prime}$, and
$\mathcal{M}\eta_{c}\eta$, respectively. To extract numerical values of
$g_{i}$, we use QCD three-point sum rule method. Predictions for the mass and
width of the molecule $\mathcal{M}$ are compared with the LHCb data for the
resonance $X(3960)$ and structure $X_{0}(4140)$. They are also confronted with
parameters of the diquark-antidiquark state
$X=[cs][\overline{c}\overline{s}]$.

This paper is organized in the following way: In Sec. II, we compute the mass
and current coupling of the molecule $\mathcal{M}$ by employing the QCD two-
point sum rule method. The processes $\mathcal{M}\rightarrow
D_{s}^{+}D_{s}^{-}$ and $\mathcal{M}\rightarrow\eta_{c}\eta^{(\prime)}$ are
considered in Sec. III, where we determine the couplings $g_{i}$ and partial
width of these decays. The full width of $\mathcal{M}$ is evaluated also in
this section. The Sec. IV is reserved for our conclusions.

##  II Spectroscopic parameters of the molecule
$\mathcal{M}=D_{s}^{+}D_{s}^{-}$

The mass and current coupling of the molecule $\mathcal{M}=D_{s}^{+}D_{s}^{-}$
can be evaluated using the QCD two-point sum rule method Shifman:1978bx ;
Shifman:1978by . This approach works quite well not only for analysis of
conventional particles, but leads to reliable predictions also in the case of
multiquark hadrons.

The principal quantity in this method is an interpolating current for a hadron
under analysis. In the case of the molecule $D_{s}^{+}D_{s}^{-}$, this current
$J(x)$ has a rather simple form

$J(x)=\left[\overline{s}_{a}(x)i\gamma_{5}c_{a}(x)\right]\left[\overline{c}_{b}(x)i\gamma_{5}s_{b}(x)\right],$
(3)

where $a$, and $b$ are color indices. This current belongs to
$\mathbf{[1}_{c}\mathbf{]}_{\overline{s}c}\mathbf{\otimes[1}_{c}\mathbf{]}_{\overline{c}s}$
representation of the color group $SU_{c}(3)$. It corresponds to a molecule
structure with spin-parities $J^{\mathrm{PC}}=0^{++}$, but may couple also to
different diquark-antidiquark structures and other four-quark hadronic
molecules Agaev:2022pis ; Chen:2022sbf ; Xin:2021wcr .

To find sum rules for the mass $m$ and coupling $f$ of the molecule
$\mathcal{M}$, one has to start from calculation of the following correlation
function

$\Pi(p)=i\int d^{4}xe^{ipx}\langle
0|\mathcal{T}\\{J(x)J^{{\dagger}}(0)\\}|0\rangle.$ (4)

At the first stage, we express $\Pi(p)$ in terms of the spectral parameters of
$\mathcal{M}$. To this end, it is necessary to insert a complete set of states
with $J^{\mathrm{PC}}=0^{++}$ into the correlation function $\Pi(p)$ and carry
out integration over $x$. These operations lead to the simple formula

$\Pi^{\mathrm{Phys}}(p)=\frac{\langle
0|J|\mathcal{M}(p\rangle\langle\mathcal{M}(p)|J^{\dagger}|0\rangle}{m^{2}-p^{2}}+\cdots.$
(5)

The expression derived by this way is a hadronic representation of the
correlator $\Pi(p)$, which forms the phenomenological side of sum rules. The
term written down explicitly in Eq. (5) is a contribution of the ground-state
particle $\mathcal{M}$, whereas contributions coming from higher resonances
and continuum states are denoted by dots.

The function $\Pi^{\mathrm{Phys}}(p)$ can be rewritten in a more convenient
form using the matrix element

$\langle 0|J|\mathcal{M}(p)\rangle=fm.$ (6)

Then, it is not difficult to find $\Pi^{\mathrm{Phys}}(p)$ in terms of the
parameters $m$ and $f$

$\Pi^{\mathrm{Phys}}(p)=\frac{m^{2}f^{2}}{m^{2}-p^{2}}+\cdots.$ (7)

The Lorentz structure of $\Pi^{\mathrm{Phys}}(p)$ has a simple form, and
consists of a term proportional to $\mathrm{I}$. Then, the invariant amplitude
$\Pi^{\mathrm{Phys}}(p^{2})$ corresponding to this structure is given by
expression in the right hand-side of Eq. (7).

The QCD side of the sum rules $\Pi^{\mathrm{OPE}}(p)$ is determined by Eq. (4)
calculated using the $c$ and $s$-quarks propagators. To this end, we insert
the explicit form of $J(x)$ into Eq. (4), contract heavy and light quark
fields, and write obtained expression using quark propagators. After these
operations, we get

$\displaystyle\Pi^{\mathrm{OPE}}(p)=i\int
d^{4}xe^{ipx}\mathrm{Tr}\left[\gamma_{5}S_{c}^{aa^{\prime}}(x)\right.$
$\displaystyle\left.\times\gamma_{5}S_{s}^{a^{\prime}a}(-x)\right]\mathrm{Tr}\left[\gamma_{5}S_{c}^{b^{\prime}b}(-x)\gamma_{5}S_{s}^{bb^{\prime}}(x)\right].$
(8)

Here $S_{c}(x)$ and $S_{s}(x)$ are the $c$ and $s$-quark propagators, explicit
expressions of which can be found in Ref. Agaev:2020zad .

The correlation function $\Pi^{\mathrm{OPE}}(p)$ is calculated by employing
the quark propagators with some fixed accuracy of the operator product
expansion ($\mathrm{OPE}$). The $\Pi^{\mathrm{OPE}}(p)$ has a trivial Lorentz
structure $\sim\mathrm{I}$ as well. Having denoted the corresponding invariant
amplitude by $\Pi^{\mathrm{OPE}}(p^{2})$ and equated it to
$\Pi^{\mathrm{Phys}}(p^{2})$, we get a sum rule equality, which can be
undergone further processing. The ground-state term and ones due to higher
resonances and continuum states contribute to this sum rule equality on equal
footing. There is a necessity to suppress unwanted contributions of higher
resonances and subtract them from this expression. For these purposes, we
apply the Borel transformation to its both sides. This operation suppresses
effects of higher resonances and continuum states, but at the same time
generates dependence of obtained equality on the Borel parameter $M^{2}$.
Afterwards, using the assumption about quark-hadron duality, we perform
continuum subtraction, which leads to additional parameter $s_{0}$ in
formulas.

The Borel transformation of the main term in $\Pi^{\mathrm{Phys}}(p^{2})$ has
a simple form

$\Pi^{\mathrm{Phys}}(M^{2})=m^{2}f^{2}e^{-m^{2}/M^{2}}.$ (9)

For the Borel transformed and subtracted amplitude
$\Pi^{\mathrm{OPE}}(p^{2})$, we find

$\Pi(M^{2},s_{0})=\int_{4(m_{c}+m_{s})^{2}}^{s_{0}}ds\rho^{\mathrm{OPE}}(s)e^{-s/M^{2}}+\Pi(M^{2}).$
(10)

The first term in Eq. (10) contains an essential part of contributions, and is
expressed using the two-point spectral density $\rho^{\mathrm{OPE}}(s)$,
derived as an imaginary part of the correlation function. The second term
$\Pi(M^{2})$ collects nonperturbative contributions extracted directly from
$\Pi^{\mathrm{OPE}}(p)$. Formulas of the functions $\rho^{\mathrm{OPE}}(s)$
and $\Pi^{\mathrm{OPE}}(p)$ are lengthy, therefore we refrain from writing
down them explicitly.

The sum rules for $m$ and $f$ read

$m^{2}=\frac{\Pi^{\prime}(M^{2},s_{0})}{\Pi(M^{2},s_{0})},$ (11)

and

$f^{2}=\frac{e^{m^{2}/M^{2}}}{m^{2}}\Pi(M^{2},s_{0}),$ (12)

where $\Pi^{\prime}(M^{2},s_{0})=d\Pi(M^{2},s_{0})/d(-1/M^{2})$.

Numerical calculations of $m$ and $f$ should be performed in accordance with
Eqs. (11) and (12), but only after fixing different vacuum condensates and
working windows for the parameters $M^{2}$ and $s_{0}$. The quark, gluon and
mixed condensates are universal and well known quantities Shifman:1978bx ;
Shifman:1978by ; Ioffe:1981kw ; Ioffe:2005ym ; Narison:2015nxh . Their
numerical values, extracted from numerous processes are listed below

$\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{s}g_{s}\sigma
Gs\rangle=m_{0}^{2}\langle\overline{s}s\rangle,\ m_{0}^{2}=(0.8\pm
0.1)~{}\mathrm{GeV}^{2},\ $
$\displaystyle\langle\frac{\alpha_{s}G^{2}}{\pi}\rangle=(0.012\pm
0.004)~{}\mathrm{GeV}^{4},$ $\displaystyle\langle
g_{s}^{3}G^{3}\rangle=(0.57\pm 0.29)~{}\mathrm{GeV}^{6},$ $\displaystyle
m_{c}=(1.27\pm 0.02)~{}\mathrm{GeV},\ \ \ m_{s}=93_{-5}^{+11}~{}\mathrm{MeV}.$
(13)

We have included the masses of the $c$ and $s$-quarks into Eq. (13) as well.
The correlation function $\Pi(M^{2},s_{0})$ is calculated by taking into
account vacuum condensates up to dimension $10$. The expression of
$\Pi(M^{2},s_{0})$ is rather lengthy, therefore we do not provide it here
explicitly. In numerical analysis we set $m_{s}^{2}=0$, but take into account
contributions proportional to $m_{s}$.

Figure 1: Dependence of the pole contribution $\mathrm{PC}$ on the Borel
parameter $M^{2}$ at different $s_{0}$. The limit $\mathrm{PC}=0.5$ is shown
by the horizontal black line. The red triangle shows the point, where the mass
$m$ of the molecule $\mathcal{M}=D_{s}^{+}D_{s}^{-}$ has effectively been
calculated.

To carry out numerical analysis one also needs to choose working regions for
the Borel and continuum subtraction parameters $M^{2}$ and $s_{0}$. They
should satisfy standard constraints of sum rule calculations. Thus, the
parameters $M^{2}$ and $s_{0}$ employed in calculations have to guarantee the
dominance of the pole contribution ($\mathrm{PC}$) and convergence of
$\mathrm{OPE}$. The former can be defined by the expression

$\mathrm{PC}=\frac{\Pi(M^{2},s_{0})}{\Pi(M^{2},\infty)},$ (14)

whereas to make sure the operator product expansion converges, we utilize the
ratio

$R(M^{2})=\frac{\Pi^{\mathrm{DimN}}(M^{2},s_{0})}{\Pi(M^{2},s_{0})},$ (15)

where $\Pi^{\mathrm{DimN}}(M^{2},s_{0})$ is a sum of a few last terms in
$\mathrm{OPE}$.

The $R(M^{2})$ and $\mathrm{PC}$ are used to restrict the lower and upper
bounds for the Borel parameter, respectively. In fact, at
$M_{\mathrm{min}}^{2}$ the function $R(M_{\mathrm{min}}^{2})$ should be less
than some fixed value, whereas at $M_{\mathrm{max}}^{2}$ the pole contribution
$\mathrm{PC}$ has to overshoot minimally acceptable limit for this parameter.
In our present analysis, we impose on $\mathrm{PC}$ and $R(M^{2})$ the
following constraints

$\mathrm{PC}\geq 0.5,\ \ R(M_{\mathrm{min}}^{2})\leq 0.05.$ (16)

The first criterion in Eq. (16) is usual for investigations of conventional
hadrons and ensures the dominance of the pole contribution. It may be used in
analysis of multiquark hadrons, as well. But in the case of multiquark
particles this constraint reduces a window for $M^{2}$. The second condition
is required to enforce the convergence of $\mathrm{OPE}$.

Dominance of a perturbative contribution to $\Pi(M^{2},s_{0})$ over a
nonperturbative term, as well as maximum stability of extracted physical
quantities when varying $M^{2}$ are among constraints to determine working
regions for the parameters $M^{2}$ and $s_{0}$. Numerical tests carried out by
taking into account these aspects of the sum rule analysis demonstrate that
the windows for $M^{2}$ and $s_{0}$,

$M^{2}\in[3,4]~{}\mathrm{GeV}^{2},\ s_{0}\in[21,22]~{}\mathrm{GeV}^{2},$ (17)

comply with aforementioned restrictions. Indeed, at
$M^{2}=4~{}\mathrm{GeV}^{2}$ the pole contribution on average in $s_{0}$
equals to $0.52$, and amounts to $0.88$ at $M^{2}=3~{}\mathrm{GeV}^{2}$. To
visualize dynamics of the pole contribution, we depict in Fig. 1 $\mathrm{PC}$
as a function of $M^{2}$ at different $s_{0}$. One can see, that the pole
contribution overshoots $0.5$ for all values of the parameters $M^{2}$ and
$s_{0}$ from Eq. (17).

To be convinced in convergence of $\mathrm{OPE}$, we calculate $\
R(M_{\mathrm{min}}^{2})$ at the minimum point
$M_{\mathrm{min}}^{2}=3~{}\mathrm{GeV}^{2}$ using in Eq. (15)
dimension-$8,9,10$ contributions. At the minimal value of $M^{2}$, we find
$R(3~{}\mathrm{GeV}^{2})\approx 0.01$ in accordance with the constraint from
Eq. (17). Results of a more detailed analysis are shown in Fig. 2. In this
figure, one sees contributions to the correlation function $\Pi(M^{2},s_{0})$
arising from the perturbative and nonperturbative terms up to dimension $8$.
The perturbative contribution forms $0.65$ part of $\Pi(M^{2},s_{0})$ at
$M^{2}=3~{}\mathrm{GeV}^{2}$ and exceeds the sum of nonperturbative terms in
the whole region of $M^{2}$. The $\mathrm{Dim3}$ term overshoots effects of
other nonperturbative operators, which enter to $\Pi(M^{2},s_{0})$ with
different signs. The $\mathrm{Dim9}$ and $\mathrm{Dim10}$ terms are
numerically very small and not demonstrated in the figure.

Our predictions for the mass $m$ and coupling $f$ read

$\displaystyle m$ $\displaystyle=$ $\displaystyle(4117\pm 85)~{}\mathrm{MeV},$
$\displaystyle f$ $\displaystyle=$ $\displaystyle(5.9\pm 0.7)\times
10^{-3}~{}\mathrm{GeV}^{4}.$ (18)

The results for $m$ and $f$ are calculated as their values averaged over the
working regions Eq. (17). Effectively they correspond to the sum rule
predictions at $M^{2}=3.5~{}\mathrm{GeV}^{2}$ and
$s_{0}=21.5~{}\mathrm{GeV}^{2}$ which is a middle point of the regions Eq.
(17): the red triangle in Fig. 1 marks exactly this point. The pole
contribution there is equal to $\mathrm{PC}\approx 0.68$, which in conjunction
with other constraints ensures ground-state nature of the molecule
$D_{s}^{+}D_{s}^{-}$ and credibility of obtained predictions.

In Fig. 3, we depict the mass $m$ of the molecule $D_{s}^{+}D_{s}^{-}$ as a
function of the Borel and continuum subtraction parameters. Physical
quantities obtained from the sum rule analysis should be stable against
variations of the Borel parameter. But the mass $m$ depends on working windows
chosen for their calculations. In fact, although the region for the Borel
parameter $M^{2}$ leads to approximately stable prediction for $m$, there is
still residual dependence on it. This effect generates theoretical
uncertainties of the sum rule calculations. It is worth noting that the
uncertainties of $m$ are smaller than ones for the coupling $f$. The reason is
that the mass $m$ is given by the ratio of the correlation functions which
compensates changes of $m$ against $M^{2}$ and $s_{0}$. The coupling $f$ , at
the same time, depends on $\Pi(M^{2},s_{0})$ and is open for an impact of the
parameters $M^{2}$ and $s_{0}$. As a result, uncertainties of calculations are
equal to $\pm 2\%$ in the case of the mass, and to $\pm 12\%$ for the current
coupling.

The region for $s_{0}$ together with $M^{2}$ has to provide the dominance of
$\mathrm{PC}$ and convergence of the operator product expansion. The parameter
$\sqrt{s_{0}}$ carries also useful information on a mass $m^{\ast}$ of the
first radial excitation of the molecule $D_{s}^{+}D_{s}^{-}$. Thus, in the
”ground-state +continuum” scheme adopted in the present work, $\sqrt{s_{0}}$
should be less than the mass $m^{\ast}$ of the first excited state. This fact
allows us to estimate low limit for $m^{\ast}\geq m+480~{}\mathrm{MeV}$.

Figure 2: Normalized contributions to $\Pi(M^{2},s_{0})$ as functions of the
Borel parameter $M^{2}$. All curves have been calculated at
$s_{0}=21.5~{}\mathrm{GeV}^{2}$.

Figure 3: Mass $m$ of the molecule $\mathcal{M}$ as a function of the Borel
$M^{2}$ (left panel), and continuum threshold $s_{0}$ parameters (right
panel).

Our result for the mass of the molecule $\mathcal{M}$ overshoots the LHCb
datum $m_{1\mathrm{exp}}$, but nicely agrees with $m_{2\mathrm{exp}}$.

##  III Full width of $\mathcal{M}=D_{s}^{+}D_{s}^{-}$

The mass and current coupling of the molecule $\mathcal{M}$ calculated in the
previous section provide information to select its possible decay modes.
Besides, one should take into account its quantum numbers
$J^{\mathrm{PC}}=0^{++}$. Because the structures $X(3960)$ and $X_{0}(4140)$
were discovered in the invariant mass distribution of the $D_{s}^{+}D_{s}^{-}$
mesons, we consider $\mathcal{M}\rightarrow D_{s}^{+}D_{s}^{-}$ as a dominant
decay mode of $\mathcal{M}$. The two-meson threshold for this decay is
$\approx 3937~{}\mathrm{MeV}$, which makes $\mathcal{M}\rightarrow
D_{s}^{+}D_{s}^{-}$ kinematically allowed channel for $\mathcal{M}$. Two other
processes which we study here are decays
$\mathcal{M}\rightarrow\eta_{c}\eta^{\prime}$ and
$\mathcal{M}\rightarrow\eta_{c}\eta$. The two-meson threshold
$3941~{}\mathrm{MeV}$ for these decays is below the mass $m$ of the molecule
$\mathcal{M}$ as well. It is easy to prove that these decays conserve the
spin, $\mathrm{P}$ and $\mathrm{C}$ parities of the initial particle
$\mathcal{M}$.

###  III.1 Process $\mathcal{M}\rightarrow D_{s}^{+}D_{s}^{-}$

The partial width of the decay $\mathcal{M}\rightarrow D_{s}^{+}D_{s}^{-}$ is
governed by a coupling $g_{1}$ that describes strong interaction at the vertex
$\mathcal{M}D_{s}^{+}D_{s}^{-}$. This partial width depends also on masses and
decay constants of the molecule $\mathcal{M}$ and mesons $D_{s}^{+}$ and
$D_{s}^{-}$. The mass and coupling of $\mathcal{M}$ have been calculated in
the present article, whereas physical parameters of the mesons $D_{s}^{+}$ and
$D_{s}^{-}$ are known from independent sources. Therefore, only physical
quantity to be found is the strong coupling $g_{1}$.

To determine $g_{1}$, we employ the QCD three-point sum rule method, and begin
our exploration from the correlation function

$\displaystyle\Pi(p,p^{\prime})=i^{2}\int
d^{4}xd^{4}ye^{i(p^{\prime}y-px)}\langle 0|\mathcal{T}\\{J^{D_{s}^{+}}(y)$
$\displaystyle\times J^{D_{s}^{-}}(0)J^{\dagger}(x)\\}|0\rangle.$ (19)

Here, $J(x)$,$\ J^{D_{s}^{+}}(y)$ and $J^{D_{s}^{-}}(0)$ are the interpolating
currents for $\mathcal{M}$, and pseudoscalar mesons $D_{s}^{+}$ and
$D_{s}^{-}$, respectively. The currents $J^{D_{s}^{+}}$ and $J^{D_{s}^{-}}$
are given by the expressions

$\displaystyle\ J^{D_{s}^{+}}(x)$ $\displaystyle=$
$\displaystyle\overline{s}_{j}(x)i\gamma_{5}c_{j}(x),$ $\displaystyle
J^{D_{s}^{-}}(x)$ $\displaystyle=$
$\displaystyle\overline{c}_{i}(x)i\gamma_{5}s_{i}(x),$ (20)

where $i$ and $j$ are color indices. The four-momenta of $\mathcal{M}$ and
$D_{s}^{+}$ are labeled by $p$ and $p^{\prime}$: Then, the momentum of the
meson $D_{s}^{-}$ is equal to $q=p-p^{\prime}$.

To find $g_{1}$, we apply standard recipes of the sum rule method and
calculate the correlation function $\Pi(p,p^{\prime})$. For these purpose, we
use the physical parameters of the molecule $\mathcal{M}$ and mesons involved
into this decay. The correlator $\Pi(p,p^{\prime})$ obtained by this manner
forms the physical side $\Pi^{\mathrm{Phys}}(p,p^{\prime})$ of the sum rule.
It is easy to see, that

$\displaystyle\Pi^{\mathrm{Phys}}(p,p^{\prime})=\frac{\langle
0|J^{D_{s}^{+}}|D_{s}^{+}(p^{\prime})\rangle\langle
0|J^{D_{s}^{-}}|D_{s}^{-}(q)\rangle}{(p^{\prime
2}-m_{D_{s}}^{2})(q^{2}-m_{D_{s}}^{2})}$ $\displaystyle\times\frac{\langle
D_{s}^{-}(q)D_{s}^{+}(p^{\prime})|\mathcal{M}(p)\rangle\langle\mathcal{M}(p)|J^{\dagger}|0\rangle}{(p^{2}-m^{2})}+\cdots,$
(21)

with $m_{D_{s}}$ being the mass of the mesons $D_{s}^{\pm}$. To get Eq. (21),
we isolate contributions of the ground-state and higher resonances and
continuum state particles from each other. In Eq. (21) the ground-state term
is written down explicitly, whereas other contributions are denoted by
ellipses.

The function $\Pi^{\mathrm{Phys}}(p,p^{\prime})$ can be rewritten in terms of
the $D_{s}^{\pm}$ mesons matrix elements

$\langle
0|J^{D_{s}^{\pm}}|D_{s}^{\pm}\rangle=\frac{m_{D_{s}}^{2}f_{D_{s}}}{m_{c}+m_{s}},\
$ (22)

where $f_{D_{s}}$ is their decay constant. We model the vertex
$\mathcal{M}D_{s}^{+}D_{s}^{-}$ by the matrix element

$\langle
D_{s}^{-}(q)D_{s}^{+}(p^{\prime})|\mathcal{M}(p)\rangle=g_{1}(q^{2})p\cdot
p^{\prime}.$ (23)

Using these expressions, it is not difficult to calculate the new expression
of the correlation function $\Pi^{\mathrm{Phys}}(p,p^{\prime})$

$\displaystyle\Pi^{\mathrm{Phys}}(p,p^{\prime})=g_{1}(q^{2})\frac{m_{D_{s}}^{4}f_{D_{s}}^{2}fm}{(m_{c}+m_{s})^{2}(p^{2}-m^{2})}$
$\displaystyle\times\frac{1}{(p^{\prime
2}-m_{D_{s}}^{2})(q^{2}-m_{D_{s}}^{2})}\frac{m^{2}+m_{D_{s}}^{2}-q^{2}}{2}+\cdots.$
(24)

The double Borel transformation of the function
$\Pi^{\mathrm{Phys}}(p,p^{\prime})$ over variables $p^{2}$ and $p^{\prime 2}$
is given by the formula

$\displaystyle\mathcal{B}\Pi^{\mathrm{Phys}}(p,p^{\prime})=g_{1}(q^{2})\frac{m_{D_{s}}^{4}f_{D_{s}}^{2}fm}{(m_{c}+m_{s})^{2}(q^{2}-m_{D_{s}}^{2})}e^{-m^{2}/M_{1}^{2}}$
$\displaystyle\times
e^{-m_{D_{s}}^{2}/M_{2}^{2}}\frac{m^{2}+m_{D_{s}}^{2}-q^{2}}{2}+\cdots.$ (25)

The correlator $\Pi^{\mathrm{Phys}}(p,p^{\prime})$ and its Borel
transformation has simple Lorentz structure $\sim\mathrm{I}$. Then the whole
expression in Eq. (24) determines the invariant amplitude
$\Pi^{\mathrm{Phys}}(p^{2},p^{\prime 2},q^{2})$.

To find the QCD side of the three-point sum rule, we express
$\Pi(p,p^{\prime})$ in terms of quark propagators, and get

$\displaystyle\Pi^{\mathrm{OPE}}(p,p^{\prime})=i^{2}\int
d^{4}xd^{4}ye^{i(p^{\prime}y-px)}$
$\displaystyle\times\mathrm{Tr}\left[\gamma_{5}S_{c}^{ia}(y-x)\gamma_{5}S_{s}^{ai}(x-y)\right]$
$\displaystyle\times\mathrm{Tr}\left[\gamma_{5}S_{s}^{jb}(-x)\gamma_{5}S_{c}^{bj}(x)\right].$
(26)

The correlator $\Pi^{\mathrm{OPE}}(p,p^{\prime})$ is computed by taking into
account nonperturbative contributions up to dimension $6$, and as
$\Pi^{\mathrm{Phys}}(p,p^{\prime})$, contains the Lorentz structure
proportional to I. Denoting relevant invariant amplitude by
$\Pi^{\mathrm{OPE}}(p^{2},p^{\prime 2},q^{2})$ , equating its double Borel
transformation $\mathcal{B}\Pi^{\mathrm{OPE}}(p^{2},p^{\prime 2},q^{2})$ to
$\mathcal{B}\Pi^{\mathrm{Phys}}(p^{2},p^{\prime 2},q^{2})$, and performing
continuum subtraction, we get the sum rule for the coupling $g_{1}(q^{2})$.

After these manipulations, the amplitude $\Pi^{\mathrm{OPE}}(p^{2},p^{\prime
2},q^{2})$ can be rewritten using the spectral density
$\rho(s,s^{\prime},q^{2})$, which is extracted as an imaginary part of
$\Pi^{\mathrm{OPE}}(p,p^{\prime})$

$\displaystyle\Pi(\mathbf{M}^{2},\mathbf{s}_{0},q^{2})=\int_{4(m_{c}+m_{s})^{2}}^{s_{0}}ds\int_{(m_{c}+m_{s})^{2}}^{s_{0}^{\prime}}ds^{\prime}\rho(s,s^{\prime},q^{2})$
$\displaystyle\times e^{-s/M_{1}^{2}}e^{-s^{\prime}/M_{2}^{2}},$ (27)

where $\mathbf{M}^{2}=(M_{1}^{2},\ M_{2}^{2})$ and $\mathbf{s}_{0}=(s_{0},\
s_{0}^{\prime})$ are the Borel and continuum threshold parameters. The couples
$(M_{1}^{2},\ s_{0})$ and $\mathbf{s}_{0}=(M_{2}^{2},\ s_{0}^{\prime})$
correspond to the molecule $M$ and $D_{s}^{+}$ meson channels, respectively.

The sum rule for $g_{1}(q^{2})$ is determined by the formula

$\displaystyle g_{1}(q^{2})$ $\displaystyle=$
$\displaystyle\frac{2(m_{c}+m_{s})^{2}}{m_{D_{s}}^{4}f_{D_{s}}^{2}fm}\frac{q^{2}-m_{D_{s}}^{2}}{m^{2}+m_{D_{s}}^{2}-q^{2}}$
(28) $\displaystyle\times
e^{m^{2}/M_{1}^{2}}e^{m_{D_{s}}^{2}/M_{2}^{2}}\Pi(\mathbf{M}^{2},\mathbf{s}_{0},q^{2}).$

Expression of $g_{1}(q^{2})$ depends on the spectroscopic parameters of the
molecule $\mathcal{M}$, and masses and decay constants of the mesons
$D_{s}^{\pm}$: they are input parameters of numerical computations. Values of
these parameters as well as masses and decay constants $f_{D_{s}}$ and
$f_{\eta_{c}}$ which are necessary to study other decays are collected in
Table 1. The masses all of mesons are borrowed from Ref. PDG:2022 . For decay
constant of the mesons $D_{s}^{\pm}$, we employ information from the same
source, whereas for $f_{\eta_{c}}$ use a prediction made in Ref.
Colangelo:1992cx  on the basis of the sum rule method.

The partial width of the decay $\mathcal{M}\rightarrow D_{s}^{+}D_{s}^{-}$
besides various input parameters, also is determined by the strong coupling
$g_{1}(m_{D_{s}}^{2})$ at the mass shell $q^{2}=m_{D_{s}}^{2}$ of the meson
$D_{s}^{-}$. At the same time, the sum rule computations of $g_{1}$ can be
carried out in deep-Euclidean region $q^{2}<0$. For simplicity, we introduce a
variable $Q^{2}=-q^{2}$ and in what follows label the obtained function by
$g_{1}(Q^{2})$. An explored range of $Q^{2}$ covers the region
$Q^{2}=1-5~{}\mathrm{GeV}^{2}$.

Numerical calculations also require choosing of the working regions for the
Borel and continuum subtraction parameters $\mathbf{M}^{2}$ and
$\mathbf{s}_{0}$. Limits imposed on $\mathbf{M}^{2}$ and $\mathbf{s}_{0}$ are
standard for sum rule calculations and were considered in the previous
section. The regions for $M_{1}^{2}$ and $s_{0}$ associated with the
$\mathcal{M}$ channel are fixed in accordance with Eq. (17). The parameters
$(M_{2}^{2},\ s_{0}^{\prime})$ for the $D_{s}^{+}$ meson channel are varied
within limits

$M_{2}^{2}\in[2.5,3.5]~{}\mathrm{GeV}^{2},\
s_{0}^{\prime}\in[5,6]~{}\mathrm{GeV}^{2}.$ (29)

Regions for $\mathbf{M}^{2}$ and $\mathbf{s}_{0}$ are chosen in such a way
that to minimize their effects of the coupling $g_{1}(Q^{2})$.

Results of computations are pictured in Fig. 4. It is seen that results for
$g_{1}(Q^{2})$ are extracted at the region $Q^{2}>0$, where sum rule gives
reliable predictions. It has just been explained above, that we need
$g_{1}(Q^{2})$ at $Q^{2}=-m_{D_{s}}^{2}$. To this end, one has to introduce
some fit function $\mathcal{G}_{1}(Q^{2})$, which at the momenta $Q^{2}>0$
gives the same values as the sum rule computations, but can be easily
extrapolated to the region $Q^{2}<0$. There are different choices for such
functions. In this article, we use $\mathcal{G}_{i}(Q^{2}),\ i=1,2,3$

$\mathcal{G}_{i}(Q^{2})=\mathcal{G}_{i}^{0}\mathrm{\exp}\left[c_{i}^{1}\frac{Q^{2}}{m^{2}}+c_{i}^{2}\left(\frac{Q^{2}}{m^{2}}\right)^{2}\right],$
(30)

where $\mathcal{G}_{i}^{0}$, $c_{i}^{1}$ and $c_{i}^{2}$ are parameters, which
should be extracted from fitting procedures. Calculations demonstrate that
$\mathcal{G}_{1}^{0}=1.29~{}\mathrm{GeV}^{-1}$, $c_{1}^{1}=2.38$, and
$c_{1}^{2}=-1.84$ lead a reasonable agreement with the sum rule’s data (see,
Fig. 4).

At the mass shell $q^{2}=m_{D_{s}}^{2}$ this function gives

$g_{1}\equiv\mathcal{G}_{1}(-m_{D_{s}}^{2})=(6.8\pm 1.6)\times 10^{-1}\
\mathrm{GeV}^{-1}.$ (31)

The width of the decay $\mathcal{M}\rightarrow D_{s}^{+}D_{s}^{-}$ is
calculated by means of the formula

$\Gamma\left[\mathcal{M}\rightarrow
D_{s}^{+}D_{s}^{-}\right]=g_{1}^{2}\frac{m_{D_{s}}^{2}\lambda}{8\pi}\left(1+\frac{\lambda^{2}}{m_{D_{s}}^{2}}\right),$
(32)

where $\lambda=\lambda\left(m,m_{D_{s}},m_{D_{s}}\right)$ and

$\displaystyle\lambda\left(a,b,c\right)$ $\displaystyle=$
$\displaystyle\frac{1}{2a}\left[a^{4}+b^{4}+c^{4}\right.$ (33)
$\displaystyle\left.-2\left(a^{2}b^{2}+a^{2}c^{2}+b^{2}c^{2}\right)\right]^{1/2}.$

Using the coupling Eq. (31), it is not difficult to compute the width of the
decay $\mathcal{M}\rightarrow D_{s}^{+}D_{s}^{-}$

$\Gamma\left[\mathcal{M}\rightarrow D_{s}^{+}D_{s}^{-}\right]=(46.5\pm
11.6)~{}\mathrm{MeV}.$ (34)

Figure 4: The sum rule results and fit function for the strong coupling $g_{1}(Q^{2})$. The red diamond denotes the point $Q^{2}=-m_{D_{s}}^{2}$.  Quantity | Value (in $\mathrm{MeV}$ units)
---|---
$m_{D_{s}}$ | $1969.0\pm 1.4$
$m_{\eta_{c}}$ | $2983.9\pm 0.4$
$m_{\eta^{\prime}}$ | $957.78\pm 0.06$
$m_{\eta}$ | $547.862\pm 0.017$
$f_{D_{s}}$ | $249.9\pm 0.5$
$f_{\eta_{c}}$ | $320\pm 40$
Table 1: Masses and decay constants of the mesons $D_{s}^{\pm}$, $\eta_{c}$,
$\eta^{\prime}$, and $\eta$ which are employed in numerical calculations.

###  III.2 Decays $\mathcal{M}\rightarrow\eta_{c}\eta^{\prime}$ and
$\mathcal{M}\rightarrow\eta_{c}\eta$

The decays $\mathcal{M}\rightarrow\eta_{c}\eta^{\prime}$ and
$\mathcal{M}\rightarrow\eta_{c}\eta$ are studied by a way described above.
But, here we take into account peculiarities of the $\eta-\eta^{\prime}$
mesons connected with mixing in this system due to $U(1)$ anomaly
Feldmann:1998vh . This effect modifies a choice of interpolating currents and
matrix elements for these particles. Although $\eta-\eta^{\prime}$ mixing can
be considered in different approaches, we use the quark-flavor basis
$|\eta_{q}\rangle=(\overline{u}u+\overline{d}d)/\sqrt{2}$ and
$|\eta_{s}\rangle=\overline{s}s$, where the physical particles $\eta$ and
$\eta^{\prime}$ have simple decompositions Feldmann:1998vh ; Agaev:2014wna ;
Agaev:2015faa

$\displaystyle\eta$ $\displaystyle=$
$\displaystyle|\eta_{q}\rangle\cos\varphi-|\eta_{s}\rangle\sin\varphi,$
$\displaystyle\eta^{\prime}$ $\displaystyle=$
$\displaystyle|\eta_{q}\rangle\sin\varphi+|\eta_{s}\rangle\cos\varphi.$ (35)

where $\varphi$ is the mixing angle in the $\\{|\eta_{q}\rangle,\
|\eta_{s}\rangle\\}$ basis. Such state-mixing implies that the same assumption
is applicable to their currents, decay constants and matrix elements.

Because in the decays $\mathcal{M}\rightarrow\eta_{c}\eta^{\prime}$,
$\eta_{c}\eta$ participate only $\overline{s}s$ components of the mesons
$\eta$ and $\eta^{\prime}$, relevant interpolating currents are given by the
expressions

$\displaystyle J^{\eta}(x)$ $\displaystyle=$
$\displaystyle-\sin\varphi\overline{s}_{j}(x)i\gamma_{5}s_{j}(x),$
$\displaystyle J^{\eta^{\prime}}(x)$ $\displaystyle=$
$\displaystyle\cos\varphi\overline{s}_{j}(x)i\gamma_{5}s_{j}(x),$ (36)

where $j$ is the color index.

We start our analysis from the process
$\mathcal{M}\rightarrow\eta_{c}\eta^{\prime}$. Then, we should consider the
correlation function

$\displaystyle\widetilde{\Pi}(p,p^{\prime})$ $\displaystyle=$ $\displaystyle
i^{2}\int d^{4}xd^{4}ye^{i(p^{\prime}y-px)}\langle
0|\mathcal{T}\\{J^{\eta_{c}}(y)$ (37) $\displaystyle\times
J^{\eta^{\prime}}(0)J^{\dagger}(x)\\}|0\rangle,$

where $\ J^{\eta_{c}}(y)$ is the interpolating current of the meson $\eta_{c}$

$\ J^{\eta_{c}}(x)=\overline{c}_{i}(x)i\gamma_{5}c_{i}(x).$ (38)

The main contribution to the correlation function
$\widetilde{\Pi}(p,p^{\prime})$ has the following form

$\displaystyle\widetilde{\Pi}^{\mathrm{Phys}}(p,p^{\prime})=\frac{\langle
0|J^{\eta_{c}}|\eta_{c}(p^{\prime})\rangle\langle
0|J^{\eta^{\prime}}|\eta^{\prime}(q)\rangle}{(p^{\prime
2}-m_{\eta_{c}}^{2})(q^{2}-m_{\eta^{\prime}}^{2})}$
$\displaystyle\times\frac{\langle\eta^{\prime}(q)\eta_{c}(p^{\prime})|\mathcal{M}(p)\rangle\langle\mathcal{M}(p)|J^{\dagger}|0\rangle}{(p^{2}-m^{2})}+\cdots,$
(39)

where the ellipses stand for contributions of higher resonances and continuum
states. The function $\widetilde{\Pi}^{\mathrm{Phys}}(p,p^{\prime})$ can be
rewritten using the the matrix elements

$\displaystyle\langle 0|J^{\eta_{c}}|\eta_{c}\rangle=\frac{m_{\eta_{c}}^{2}\
f_{\eta_{c}}}{2m_{c}},\ $ $\displaystyle
2m_{s}\langle\eta^{\prime}|\overline{s}i\gamma_{5}s|0\rangle=h_{\eta^{\prime}}^{s},$
(40)

where $m_{\eta_{c}}$ and $\ f_{\eta_{c}}$ are the mass and decay constant of
the $\eta_{c}$ meson. The matrix element of local operator
$\overline{s}i\gamma_{5}s$ sandwiched between the meson $\eta^{\prime}$ and
vacuum states is denoted by $h_{\eta^{\prime}}^{s}$ Agaev:2014wna . The
parameter $h_{\eta^{\prime}}^{s}$ follows the $\eta-\eta^{\prime}$ state-
mixing pattern, and we get

$h_{\eta^{\prime}}^{s}=h_{s}\cos\varphi.$ (41)

The parameter $h_{s}$ can be defined theoretically Agaev:2014wna , but for our
purposes it is enough to use values of $h_{s}$ and $\varphi$ extracted from
phenomenological analyses

$\displaystyle h_{s}$ $\displaystyle=$ $\displaystyle(0.087\pm
0.006)~{}\mathrm{GeV}^{3},$ $\displaystyle\varphi$ $\displaystyle=$
$\displaystyle 39.3^{\circ}\pm 1.0^{\circ}.$ (42)

The vertex $M\eta_{c}\eta^{\prime}$ has the following form

$\langle\eta^{\prime}(q)\eta_{c}(p^{\prime})|M(p)\rangle=g_{2}(q^{2})p\cdot
p^{\prime},$ (43)

where $g_{2}$ is the strong coupling at the vertex $M\eta_{c}\eta^{\prime}$.
By employing these matrix elements, we obtain a new expression for
$\widetilde{\Pi}^{\mathrm{Phys}}(p,p^{\prime})$

$\displaystyle\widetilde{\Pi}^{\mathrm{Phys}}(p,p^{\prime})=g_{2}(q^{2})\frac{fmm_{\eta_{c}}^{2}\
f_{\eta_{c}}h_{s}\cos^{2}\varphi}{4m_{c}m_{s}{}(p^{2}-m^{2})}$
$\displaystyle\times\frac{1}{(p^{\prime
2}-m_{\eta_{c}}^{2})(q^{2}-m_{\eta^{\prime}}^{2})}\frac{m^{2}+m_{\eta_{c}}^{2}-q^{2}}{2}+\cdots.$

The QCD side of the sum rule for $g_{2}(q^{2})$ is given by the formula

$\displaystyle\widetilde{\Pi}^{\mathrm{OPE}}(p,p^{\prime})=-\cos\varphi\int
d^{4}xd^{4}ye^{i(p^{\prime}y-px)}$
$\displaystyle\times\mathrm{Tr}\left[\gamma_{5}S_{c}^{ia}(y-x)\gamma_{5}S_{s}^{ai}(x-y)\right]$
$\displaystyle\times\mathrm{Tr}\left[\gamma_{5}S_{s}^{jb}(-x)\gamma_{5}S_{c}^{bj}(x-y)\right],$
(45)

The sum rule for the coupling $g_{2}(q^{2})$ is obtained using Borel
transformations of invariant amplitudes
$\widetilde{\Pi}^{\mathrm{Phys}}(p^{2},p^{\prime 2},q^{2})$ and
$\widetilde{\Pi}^{\mathrm{OPE}}(p^{2},p^{\prime 2},q^{2})$ and is equal to

$\displaystyle g_{2}(q^{2})$ $\displaystyle=$
$\displaystyle-\frac{8m_{c}m_{s}}{fmm_{\eta_{c}}^{2}\
f_{\eta_{c}}h_{s}\cos\varphi}\frac{q^{2}-m_{\eta_{c}}^{2}}{m^{2}+m_{\eta_{c}}^{2}-q^{2}}$
(46) $\displaystyle\times
e^{m^{2}/M_{1}^{2}}e^{m_{\eta_{c}}^{2}/M_{2}^{2}}\widetilde{\Pi}(\mathbf{M}^{2},\mathbf{s}_{0},q^{2}).$

Here, $\widetilde{\Pi}(\mathbf{M}^{2},\mathbf{s}_{0},q^{2})$ is the Borel
transformed and subtracted amplitude
$\widetilde{\Pi}^{\mathrm{OPE}}(p^{2},p^{\prime 2},q^{2})$.

The coupling $g_{1}(q^{2})$ is calculated using the following Borel and
continuum threshold parameters in the $\eta_{c}$ channel

$M_{2}^{2}\in[3,4]~{}\mathrm{GeV}^{2},\
s_{0}^{\prime}\in[9.5,10.5]~{}\mathrm{GeV}^{2},$ (47)

whereas as $M_{1}^{2}$ and $s_{0}$ for the $M$ channel, we use Eq. (17). The
strong coupling $g_{2}$ is defined at the mass shell of the $\eta^{\prime}$
meson. The fit function $\mathcal{G}_{2}(Q^{2})$ given by Eq. (30) has the
parameters $\mathcal{G}_{2}^{0}=0.21~{}\mathrm{GeV}^{-1}$, $c_{2}^{1}=5.08$,
and $c_{2}^{2}=-4.04$. Computations yield

$g_{2}\equiv\mathcal{G}_{2}(-m_{\eta^{\prime}}^{2})=(1.6\pm 0.3)\times
10^{-1}\ \mathrm{GeV}^{-1}.$ (48)

The partial width of this decay can be evaluated using the formula Eq. (32),
in which one should make substitutions $g_{1}\rightarrow g_{2}$,
$m_{D_{s}}^{2}\rightarrow m_{\eta_{c}}^{2}$ and
$\lambda\left(m,m_{D_{s}},m_{D_{s}}\right)\rightarrow\widetilde{\lambda}\left(m,m_{\eta_{c}},m_{\eta^{\prime}}\right)$.
Then, for the width of the decay $M\rightarrow\eta_{c}\eta^{\prime}$, we find

$\Gamma\left[\mathcal{M}\rightarrow\eta_{c}\eta^{\prime}\right]=(4.9\pm
1.1)~{}\mathrm{MeV}.$ (49)

Analysis of the decay $\mathcal{M}\rightarrow\eta_{c}\eta$ can be performed in
a similar manner. Avoiding further details, let us write down predictions
obtained for key quantities. Thus, the strong coupling $g_{3}$ at the vertex
$\mathcal{M}\eta_{c}\eta$ is given by the formula

$g_{3}\equiv|\mathcal{G}_{3}(-m_{\eta}^{2})|=(1.5\pm 0.3)\times 10^{-1}\
\mathrm{GeV}^{-1},$ (50)

where parameters of the fit function are
$\mathcal{G}_{3}^{0}=-0.17~{}\mathrm{GeV}^{-1}$, $c_{2}^{1}=7.29$, and
$c_{2}^{2}=-7.33$. The width of the decay $\mathcal{M}\rightarrow\eta_{c}\eta$
is

$\Gamma\left[\mathcal{M}\rightarrow\eta_{c}\eta\right]=(7.7\pm
1.8)~{}\mathrm{MeV}.$ (51)

Then, it is not difficult to find the full width of the molecule $\mathcal{M}$

$\Gamma_{\mathcal{M}}=(59\pm 12)~{}\mathrm{MeV}.$ (52)

The width of the hadronic molecule $\mathcal{M}$, within errors of
calculations and measurements, agrees with the LHCb datum from Eq. (2).

##  IV Conclusions

In this work, we have calculated the mass and width the scalar molecule
$\mathcal{M}=D_{s}^{+}D_{s}^{-}$ in the framework of the QCD sum rule methods.
The mass of $\mathcal{M}$ has been evaluated using the two-point sum rule
approach. The full width of $\mathcal{M}$ has been computed by taking into
account decay modes $\mathcal{M}\rightarrow D_{s}^{+}D_{s}^{-}$, and
$\mathcal{M}\rightarrow\eta_{c}\eta^{(\prime)}$. Strong couplings $g_{i}$ that
determine width of these decays have been found in the framework of the three-
point sum rule method.

Our result for the mass $m=(4117\pm 85)~{}\mathrm{MeV}$ of the molecule
$\mathcal{M}$ exceeds considerably the corresponding LHCb datum
$m_{1\mathrm{exp}}$, but is comparable with $m_{2\mathrm{exp}}$. The full
width $\Gamma_{\mathcal{M}}=(59\pm 12)~{}\mathrm{MeV}$ of $\mathcal{M}$ is
consistent with the LHCb measurement $\Gamma_{2\mathrm{exp}}$ as well. It is
evident, that the molecule $\mathcal{M}$ is considerably heavier than the
resonance $X(3960)$, which makes problematic its interpretation as $X(3960)$.

The resonance $X(3960)$ was examined in our article Agaev:2022pis  as the
tetraquark $X=[cs][\overline{c}\overline{s}]$ with quantum numbers
$J^{\mathrm{PC}}=0^{++}$. We obtained the following predictions for the
parameters of $X$: the mass $m=(3976\pm 85)~{}\mathrm{MeV}$ and width
$\Gamma_{\mathrm{X}}=(42.2\pm 8.3)~{}\mathrm{MeV}$. The parameters of the
diquark-antidiquark state $X$ are in nice agreements with the LHCb data given
by Eq. (1), therefore in Ref. Agaev:2022pis , $X$ was identified with
$X(3960)$.

In the context of the sum rule approach the molecule $D_{s}^{+}D_{s}^{-}$ was
also studied in Ref. Xin:2022bzt . In accordance with this paper, the mass of
such hadronic molecule is equal to $(3980\pm 100)~{}\mathrm{MeV}$ and agrees
with the LHCb data. It is worth noting that the authors did not analyze
quantitatively the width of this state. Our results for parameters of
$\mathcal{M}$, even within existing errors of calculations, does not support
molecule assignment for $X(3960)$.

By taking into accout predictions for the mass $m$ and full width
$\Gamma_{\mathcal{M}}$ of the molecule $\mathcal{M}=D_{s}^{+}D_{s}^{-}$
obtained in the present work, we argue that the molecule $\mathcal{M}$ may be
a candidate to the structure $X_{0}(4140)$.

## References

  * (1) R. Aaij et al. [LHCb Collaboration], Phys. Rev. Lett. 118, 022003 (2017); R. Aaij et al. [LHCb Collaboration], Phys. Rev. D 95, 012002 (2017).
  * (2) T. Aaltonen et al. [CDF Collaboration], Phys. Rev. Lett. 102, 242002 (2009).
  * (3) S. Chatrchyan et al. [CMS Collaboration], Phys. Lett. B 734, 261 (2014).
  * (4) V. M. Abazov et al. [D0 Collaboration], Phys. Rev. D 89, 012004 (2014).
  * (5) R. Aaij et al. [LHCb Collaboration], Phys. Rev. Lett. 127, 082001 (2021). LHCb:2022vsv
  * (6) The LHCb Collaboration, arXiv:2210.15153 [hep-ex].
  * (7) J. Nieves, and M. P. Valderrama, Phys. Rev. D 86, 056004 (2012).
  * (8) Z. G. Wang, Eur. Phys. J. C 74, 2874 (2014).
  * (9) R. F. Lebed, and A. D. Polosa, Phys. Rev. D 93, 094024 (2016).
  * (10) W. Chen, H. X. Chen, X. Liu, T. G. Steele, and S. L. Zhu, Phys. Rev. D 96, 114017 (2017).
  * (11) L. Meng, B. Wang, and S. L. Zhu, Sci. Bull. 66, 1288 (2021).
  * (12) S. S. Agaev, K. Azizi and H. Sundu, Phys. Rev. D 95, 114003 (2017).
  * (13) H. Sundu, S. S. Agaev and K. Azizi, Phys. Rev. D 98, 054021 (2018).
  * (14) S. S. Agaev, K. Azizi and H. Sundu, Phys. Rev. D 106, 014025 (2022).
  * (15) M. Bayar, A. Feijoo, and E. Oset, arXiv:2207.08490 [hep-ph].
  * (16) T. Ji, X. K. Dong, M. Albaladejo, M. L. Du, F. K. Guo, and J.Nieves, arXiv:2207.08563 [hep-ph].
  * (17) Q. Xin, Z. G. Wang, and X. S. Yang, arXiv:2207.09910 [hep-ph].
  * (18) J. M. Xie, M. Z. Liu, and L. S. Geng, arXiv:2207.12178 [hep-ph].
  * (19) R. Chen, and Q.Huang, arXiv:2209.05180 [hep-ph].
  * (20) D. Guo, J. Z. Wang, D. Y. Chen, and X. Liu, arXiv:2210.16720 [hep-ph].
  * (21) S. S. Agaev, K. Azizi and H. Sundu, arXiv:2211.14129 [hep-ph].
  * (22) M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B 147, 385 (1979).
  * (23) M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B 147, 448 (1979).
  * (24) H. X. Chen, Y. X. Yan, W. Chen, Phys. Rev. D 106, 094019 (2022).
  * (25) Q. Xin, and Z. G. Wang, Eur. Phys. J. A 58, 118 (2022).
  * (26) S. S. Agaev, K. Azizi and H. Sundu, Turk. J. Phys. 44, 95 (2020).
  * (27) B. L. Ioffe, Nucl. Phys. B 188, 317 (1981) [erratum: Nucl. Phys. B 191, 591 (1981) ].
  * (28) B. L. Ioffe, Prog. Part. Nucl. Phys. 56, 232 (2006).
  * (29) S. Narison, Nucl. Part. Phys. Proc. 270-272, 143 (2016).
  * (30) R. L. Workman et al. [Particle Data Group], Prog. Theor. Exp. Phys. 2022, 083C01 (2022).
  * (31) P. Colangelo, G. Nardulli and N. Paver, Z. Phys. C 57, 43 (1993).
  * (32) T. Feldmann, P. Kroll, and B. Stech Phys. Rev. D 58, 114006 (1998).
  * (33) S. S. Agaev, V. M. Braun, N. Offen, F. A. Porkert and A. Schäfer, Phys. Rev. D 90, 074019 (2014).
  * (34) S. S. Agaev, K. Azizi and H. Sundu, Phys. Rev. D 92, 116010 (2015).