# Modeling the resonance $T_{cs0}^{a}(2900)^{++}$ as a hadronic molecule
$D^{\ast+}K^{\ast+}$

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, Türkiye  H. Sundu  Department of
Physics, Kocaeli University, 41380 Izmit, Türkiye  Department of Physics
Engineering, Istanbul Medeniyet University, 34700 Istanbul, Türkiye

###### Abstract

The doubly charged scalar resonance $T_{cs0}^{a}(2900)^{++}$ is studied in the
context of the hadronic molecule model. We consider $T_{cs0}^{a}(2900)^{++}$
as a molecule $M=D^{\ast+}K^{\ast+}$ composed of vector mesons, and calculate
its mass, current coupling and full width. The spectroscopic parameters of
$M$, i.e., its mass and current coupling, are found by means of the QCD two-
point sum rule method by taking into account vacuum expectation values of
quark, gluon and mixed operators up to dimension $10$. The width of the
molecule $M$ is evaluated through the calculations of the partial widths of
the decay channels $M\to D_{s}^{+}\pi^{+}$, $M\to D_{s}^{\ast+}\rho^{+}$, and
$M\to D^{\ast+}K^{\ast+}$. Partial widths of these processes are determined by
strong couplings $g_{1}$, $g_{2}$, and $g_{3}$ of particles at vertices
$MD_{s}^{+}\pi^{+}$, $MD_{s}^{\ast+}\rho^{+}$, and $MD^{\ast+}K^{\ast+}$,
respectively. We calculate the couplings $g_{i}$ by employing the QCD light-
cone sum rule approach and technical tools of the soft-meson approximation.
Predictions obtained for the mass $m=(2924\pm 107)~{}\mathrm{MeV}$ and width
$\Gamma=(123\pm 25)~{}\mathrm{MeV}$ of the hadronic molecule $M$ allow us to
consider it as a possible candidate of the resonance $T_{cs0}^{a}(2900)^{++}$.

##  I Introduction

Recently, the LHCb collaboration discovered new resonances
$T_{cs0}^{a}(2900)^{0/++}$ (in what follows, $T_{cs0}^{a0/++}$) in the
processes $B^{0}\rightarrow\overline{D}^{0}D_{s}^{+}\pi^{-}$ and
$B^{+}\rightarrow D^{-}D_{s}^{+}\pi^{+}$ LHCb:2022xob ; LHCb:2022bkt ,
respectively. They were fixed in the $D_{s}^{+}\pi^{-}$ and $D_{s}^{+}\pi^{+}$
mass distributions, and are structures with spin-parity
$J^{\mathrm{P}}=0^{+}$. From the analysis of the decay channels of
$T_{cs0}^{a0/++}$ it becomes clear that they are fully open flavor four-quark
systems of $cd\overline{s}\overline{u}$/$cu\overline{s}\overline{d}$. Their
resonant parameters are consistent with each other, which means that they are
members of an isospin triplet: This is the first observation of an isospin
triplet of exotic mesons with four different quark flavors. The resonance
$T_{cs0}^{a++}$ has an additional attractive feature as the first doubly
charged exotic meson discovered experimentally.

It should be emphasized that $T_{cs0}^{a0/++}$ are not first fully open flavor
resonances seen by the LHCb experiment. Indeed, previously LHCb informed about
scalar $X_{0}(2900)$ and vector $X_{1}(2900)$ structures (hereafter $X_{0}$
and $X_{1}$, respectively), which were found in the $D^{-}K^{+}$ invariant
mass distribution of the decay $B^{+}\rightarrow D^{+}D^{-}K^{+}$ LHCb:2020A ;
LHCb:2020 . In a four-quark picture both $X_{0}$ and $X_{1}$ have the same
contents $ud\overline{s}\overline{c}$. New resonances $T_{cs0}^{a0/++}$ fill
up the list of such particles.

The exotic mesons built of four different quarks have already attracted the
interest of researches. Relevant activities started from announcement by the
$D0$ collaboration about the resonance $X(5568)$ D0:2016mwd ; D0:2017qqm
presumably composed of quarks $su\overline{b}\overline{d}$. Despite the fact
that LHCb, CMS and ATLAS experiments did not confirm existence of this state,
technical tools elaborated during this activity led to some interesting
results, and are still in use in numerous research works. One of such results
is prediction of a charmed partner $X_{c}=[su][\overline{c}\overline{d}]$ of
$X(5568)$ in the diquark-antidiquark model Agaev:2016lkl ; Chen:2016mqt . In
our paper Agaev:2016lkl  it was investigated in a rather detailed form. Thus,
we calculated the mass and full width of this tetraquark in the context of QCD
sum rule method using different interpolating currents. In the case of scalar-
scalar current we obtained $m_{\mathrm{S}}=(2634\pm 62)~{}\mathrm{MeV}$ and$\
\Gamma_{\mathrm{S}}=(57.7\pm 11.6)~{}\mathrm{MeV}$, whereas the axial-axial
current led to predictions $m_{\mathrm{A}}=(2590\pm 60)~{}\mathrm{MeV}$ and
$\Gamma_{\mathrm{A}}=(63.4\pm 14.2)~{}\mathrm{MeV}$. It is worth noting that
an estimation $(2.55\pm 0.09)~{}\mathrm{GeV}$ for the mass of $X_{c}$ was made
in Ref. Chen:2016mqt , as well.

The doubly charged exotic mesons were also objects of interesting analyses.
The $-2|e|$ charged scalar, pseudoscalar and axial-vector diquark-antidiquarks
$Z_{\overline{c}s}=[sd][\overline{u}\overline{c}]$ were explored in Ref.
Agaev:2017oay . Another class of tetraquarks
$Z^{++}=[cu][\overline{s}\overline{d}]$ with the electric charge $2|e|$ are
antiparticles of the states $Z_{\overline{c}s}$ and have the same masses and
decay widths. Parameters of the vector tetraquark $Z_{\mathrm{V}}^{++}$ from
this group of particles were found in Ref. Agaev:2021jsz .

The discovery of the resonances $X_{0}$ and $X_{1}$ highly intensified
investigations of fully open flavor structures Karliner:2020vsi ; Wang:2020xyc
; Chen:2020aos ; Liu:2020nil ; Liu:2020orv ; Molina:2020hde ; Hu:2020mxp ;
He:2020jna ; Lu:2020qmp ; Zhang:2020oze ; Huang:2020ptc ; Xue:2020vtq ;
Yang:2021izl ; Wu:2020job ; Abreu:2020ony ; Wang:2020prk ; Xiao:2020ltm ;
Dong:2020rgs ; Burns:2020xne ; Bondar:2020eoa ; Chen:2020eyu ;
Albuquerque:2020ugi ; Wang:2021lwy . In these articles different models were
suggested to explain the observed parameters of these states and understand
their inner organizations. Traditionally they were explored in the diquark-
antidiquark and hadronic molecule pictures, which are dominant models to
account for similar experimental data. Thus, $X_{0}$ was treated as a scalar
diquark-antidiquark state $[sc][\overline{u}\overline{d}]$ in Refs.
Karliner:2020vsi ; Wang:2020xyc . The $X_{0}$ was assigned to be the $S$-wave
hadronic molecule $D^{\ast-}K^{\ast+}$, whereas $X_{1}$ was examined as the
$P$-wave diquark-antidiquark state $[ud][\overline{c}\overline{s}]$ in Ref.
Chen:2020aos . There were attempts to consider these structures as
rescattering effects. In fact, in Ref. Liu:2020orv  it was asserted that two
resonance-like peaks in the process $B^{+}\rightarrow D^{+}D^{-}K^{+}$ may be
generated by rescattering effects and occur in the LHCb experiment as the
states $X_{0}$ and $X_{1}$.

The structures $X_{0}$ and $X_{1}$ were studied in our publications as well.
The mass and width of the resonance $X_{0}$ were calculated in Ref.
Agaev:2020nrc  in the framework of a hadronic molecule model
$\overline{D}^{\ast 0}K^{\ast 0}$. Results found in this work for the
parameters of $X_{0}$ allowed us to confirm its molecule nature. We explored
also the resonance $X_{1}$ by considering it as a vector diquark-antidiquark
state $X_{\mathrm{V}}=[ud][\overline{c}\overline{s}]$ Agaev:2021knl : It
turned out that the diquark-antidiquark structure is an appropriate model to
explain the measured parameters of $X_{1}$.

The LHCb observed only the vector tetraquark
$X_{\mathrm{V}}=[ud][\overline{c}\overline{s}]$, which was interpreted as
$X_{1}$. It is quite possible that, in near future, the diquark-antidiquark
structures $[ud][\overline{c}\overline{s}]$ with other quantum numbers will be
seen in various exclusive processes. Therefore, parameters of these yet
hypothetical exotic mesons are necessary to form a theoretical basis for
upcoming experimental activities. Motivated by this reason, we computed the
masses and full widths of the ground-state and radially excited scalar
particles $X_{0}^{(\prime)}=[ud][\overline{c}\overline{s}]$ Agaev:2022eeh .
The axial-vector and pseudoscalar tetraquarks $X_{\mathrm{AV}}$ and
$X_{\mathrm{PS}}$ were investigated in Ref. Sundu:2022kyd , in which we
evaluated their spectroscopic parameters, i.e., their masses and current
couplings, and widths.

The resonances $T_{cs0}^{a0/++}$ are last experimentally confirmed members of
fully open flavor tetraquark family. In this article we are going to study the
doubly charged state $T_{cs0}^{a++}$, therefore, below, write down its
parameters measured by LHCb LHCb:2022bkt :

$\displaystyle M_{\mathrm{exp}}$ $\displaystyle=$ $\displaystyle(2921\pm 17\pm
20)~{}\mathrm{MeV},$ $\displaystyle\Gamma_{\mathrm{exp}}$ $\displaystyle=$
$\displaystyle(137\pm 32\pm 17)~{}\mathrm{MeV}.$ (1)

Observation of new tetraquarks $T_{cs0}^{a0/++}$ generated theoretical
investigations aimed to bring them under one of existing models of four-quark
mesons. In our article Agaev:2022duz , we argued that diquark-antidiquark
structures are not suitable for these resonances, because parameters of such
states were already evaluated and predictions obtained for their masses are
well below the LHCb data. One of possible ways to explain $T_{cs0}^{a0/++}$ is
to treat them as hadronic molecules. Then $T_{cs0}^{a++}$ may be interpreted
as a hadronic molecule $D_{s}^{\ast+}\rho^{+}$ or $D^{\ast+}K^{\ast+}$. In
Ref. Agaev:2022duz , we realized first of these scenarios, and estimated the
mass of molecule $D_{s}^{\ast+}\rho^{+}$ by employing the QCD two-point sum
rule approach. Our result $m=(2917\pm 135)~{}\mathrm{MeV}$ for the mass of the
molecule $D_{s}^{\ast+}\rho^{+}$ is consistent with Eq. (1).

The resonances $T_{cs0}^{a0/++}$ were investigated in Refs. Chen:2022svh ;
Ge:2022dsp ; Wei:2022wtr ; Molina:2022jcd ; Liu:2022hbk ; Yue:2022mnf ;
Qin:2022nof  as well, in which authors used different models and calculational
schemes. The one-boson exchange model was employed in Ref. Chen:2022svh  to
study interactions in systems of $D^{(\ast)}K^{(\ast)}$ mesons. Analysis
allowed the authors to assign $T_{cs0}^{a++}$ to be an isovector
$D^{\ast+}K^{\ast+}$ molecule state with the spin-parity
$J^{\mathrm{P}}=0^{+}$ and mass $2891~{}\mathrm{MeV}$. Interpretation of new
tetraquark candidate $T_{cs0}^{a}$ as the resonance-like structure generated
by threshold effects was proposed in Ref. Ge:2022dsp . It was argued that the
triangle singularity induced by the $\chi_{c1}K^{\ast}D^{\ast}$ loop peaks
around the threshold $D^{\ast}K^{\ast}$ and may simulate $T_{cs0}^{a}$. A
multiquark color flux-tube model was used to investigate the resonances
$T_{cs0}^{a0/++}$ in the context of the diquark-antidiquark model Wei:2022wtr
. The authors found that a system $[cu][\overline{s}\overline{d}]$ built of
color antitriplet diquark and triplet antidiquark with the mass
$2923~{}\mathrm{MeV}$ is a very nice candidate to the resonance
$T_{cs0}^{a++}$. Features of the charmed-strange tetraquarks were explored
also in Ref. Liu:2022hbk  in a nonrelativistic potential quark model. Decays
of the neutral state $T_{cs0}^{a0}$ in the molecular picture were considered
in Ref. Yue:2022mnf , whereas production mechanisms of the hidden- and open-
charm tetraquarks in $B$ decays were addressed in Ref. Qin:2022nof .

In the present work, we explore the resonance $T_{cs0}^{a++}$ in the context
of the hadronic molecule model. We implement the second scenario and model
$T_{cs0}^{a++}$ as the hadronic molecule $M=D^{\ast+}K^{\ast+}$. Our analyses
of $M$ include calculations of its mass $m$, current coupling $f$ and full
width $\Gamma$. The spectroscopic parameters of $M$ are extracted from the QCD
two-point sum rule computations Shifman:1978bx ; Shifman:1978by  by taking
into account vacuum expectation values of different quark, gluon and mixed
operators up to dimension $10$.

To estimate the full width of the molecule $M$, we consider its decays to
pairs of conventional mesons $D_{s}^{+}\pi^{+}$, $D_{s}^{\ast+}\rho^{+}$, and
$D^{\ast+}K^{\ast+}$. Partial widths of these processes depend on parameters
of initial and final-state particles, as well as on couplings $g_{i}$, which
determine strong interactions of the molecule $M$ and mesons at the vertices
$MD_{s}^{+}\pi^{+}$, $MD_{s}^{\ast+}\rho^{+}$, and $MD^{\ast+}K^{\ast+}$,
respectively. Because masses and decay constants of ordinary mesons are known,
and spectroscopic parameters of $M$ are object of present studies, only missed
quantities are strong couplings $g_{i}$. The couplings $g_{i}$ are evaluated
in the framework of QCD light-cone sum rule (LCSR) method Balitsky:1989ry .
Additionally, to treat technical peculiarities of four-quark molecule-two
conventional meson vertices, we invoke the technical methods of soft-meson
approximation Ioffe:1983ju ; Belyaev:1994zk .

This article is organized in the following way: In Sec. II, we find the sum
rules for the mass $m$ and current coupling $f$ of the molecule
$M=D^{\ast+}K^{\ast+}$ in the framework of QCD sum rule method. Numerical
analysis of the quantities $m$ and $f$ is carried out in this section as well,
where their values are evaluated. In section III, we investigate the vertices
$MD_{s}^{+}\pi^{+}$, $MD_{s}^{\ast+}\rho^{+}$, and $MD^{\ast+}K^{\ast+}$ and
calculate the corresponding strong couplings $g_{i},\ i=1,2,3$. Obtained
information on $g_{i}$ is used to find partial widths of these decay channels,
and estimate full width of the molecule $M$. The section IV is reserved for
our conclusions.

##  II Mass $m$ and current coupling $f$ of the hadronic molecule
$M=D^{\ast+}K^{\ast+}$

The key quantity necessary to investigate the spectroscopic parameters of the
molecule $M$ using the QCD two-point sum rule method, is the interpolating
current $J(x)$ for this state. An analytic form of $J(x)$ depends on the
structure and constituents of a four-quark exotic meson
$\overline{q}_{1}\overline{q}_{2}q_{3}q_{4}$. In the molecule picture the
color-singlet structures come from
$\mathbf{[1}_{c}\mathbf{]}_{\overline{q}_{1}q_{3}}\mathbf{\otimes[1}_{c}\mathbf{]}_{\overline{q}_{2}q_{4}}$
and
$\mathbf{[8}_{c}\mathbf{]}_{\overline{q}_{1}q_{3}}\mathbf{\otimes[8}_{c}\mathbf{]}_{\overline{q}_{2}q_{4}}$
and ($q_{3}\leftrightarrow q_{4}$) terms of the color group $SU_{c}(3)$. In
the case under discussion, we assume that the hadronic molecule $M$ is
composed of two ordinary vector mesons $D^{\ast+}$ and $K^{\ast+}$, and
restrict our analysis by singlet-singlet type current. Then, in the
$\mathbf{[1}_{c}\mathbf{]}_{\overline{d}c}\mathbf{\otimes[1}_{c}\mathbf{]}_{\overline{s}u}$
representation $J(x)$ takes the following form

$J(x)=[\overline{d}_{a}(x)\gamma^{\mu}c_{a}(x)][\overline{s}_{b}(x)\gamma_{\mu}u_{b}(x)],$
(2)

where $a$ and $b$ are color indices.

This current has the meson-meson structure and is a local product of two
vector currents corresponding to the mesons $D^{\ast+}$ and $K^{\ast+}$. It
couples well to the molecule state $D^{\ast+}K^{\ast+}$. But, at the same
time, $J(x)$ couples also to diquark-antidiquark states, because using Fierz
transformation a molecule current can be presented as a weighted sum of
different diquark-antidiquark currents Wang:2020rcx . In its turn, a diquark-
antidiquark current is expressible via molecule structures (for instance, see
Refs. Chen:2022sbf ; Xin:2021wcr ). For example, the current $J(x)$ rewritten
in the form

$J(x)=\delta_{am}\delta_{bn}[\overline{d}_{a}\gamma^{\mu}c_{m}][\overline{s}_{b}\gamma_{\mu}u_{n}],$
(3)

after Fierz transformation contains the vector-vector component

$J_{\mathrm{VV}}=-\frac{1}{2}\delta_{am}\delta_{bn}[\overline{d}_{a}\gamma^{\mu}u_{n}][\overline{s}_{b}\gamma_{\mu}c_{m}].$
(4)

Using the rearrangement of the color indices
$\delta_{am}\delta_{bn}=\delta_{an}\delta_{bm}+\epsilon_{abk}\epsilon_{mnk},$
with $\epsilon_{ijk}$ being the Levi-Civita epsilon, it is not difficult to
see that

$J_{\mathrm{VV}}=-\frac{1}{2}[\overline{d}_{a}\gamma^{\mu}u_{a}][\overline{s}_{b}\gamma_{\mu}c_{b}]+\cdots.$
(5)

The term in Eq. (5) is the $-D_{s}^{\ast+}\rho^{+}/2$ meson-meson current,
whereas dots indicate the second component which should be further manipulated
to become diquark type current(s). In other words, $J(x)$ couples also to a
molecule $D_{s}^{\ast+}\rho^{+}$ and can couple to other meson pairs such as
$D_{s}^{+}\pi^{+}$ with the same contents and quantum numbers. Nevertheless,
the current $J(x)$ corresponds mainly to the molecule $D^{\ast+}K^{\ast+}$,
which will be demonstrated quantitatively in the next section by comparing its
strong couplings $g_{i}$ to different two-meson states.

The sum rules for the mass $m$ and current coupling $f$ of the hadronic
molecule $M$ can be obtained from analysis of the correlation function
Shifman:1978bx ; Shifman:1978by ,

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

with $\mathcal{T}$ being the time-ordering operator.

To derive the required sum rules, the correlator $\Pi(p)$ should be expressed
using the physical parameters of the molecule $M$, as well as calculated in
terms of the fundamental parameters of QCD in quark-gluon language. The first
expression establishes the physical (phenomenological) side of the sum rules,
for which we get

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

To obtain $\Pi^{\mathrm{Phys}}(p)$, we insert a complete set of the
intermediate states with the content and quantum numbers of the state $M$ into
Eq. (6), and carry out integration over $x$. In Eq. (7), the contribution of
the ground-state particle $M$ is isolated and shown explicitly, whereas dots
denote effects due to higher resonances and continuum states in the $M$
channel.

For further simplification of $\Pi^{\mathrm{Phys}}(p)$, it is convenient to
introduce the physical parameters of $M$ by means of the matrix element

$\langle 0|J|M\rangle=fm.$ (8)

Then, we get the final formula for the function $\Pi^{\mathrm{Phys}}(p)$:

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

The r.h.s. of Eq. (9) contains only a trivial Lorentz structure, which is the
unit matrix $\mathrm{I}$. The function $f^{2}m^{2}/(m^{2}-p^{2})$ is the
invariant amplitude $\Pi^{\mathrm{Phys}}(p^{2})$ that corresponds to this
structure: It will be used in our following analysis.

The QCD side of the sum rules, $\Pi^{\mathrm{OPE}}(p)$, has to be calculated
in the operator product expansion ($\mathrm{OPE}$) with some fixed accuracy.
To derive $\Pi^{\mathrm{OPE}}(p)$, we insert the interpolating current $J(x)$
into Eq. (6), contract the corresponding heavy and light quark fields, and
write the obtained expression in terms of the corresponding quark propagators.
Having carried out these manipulations, we get for $\Pi^{\mathrm{OPE}}(p)$,

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

where $S_{c}(x)$ and $S_{u(s,d)}(x)$ are the quark propagators. The explicit
expressions of the heavy and light quarks propagators can be found in Ref.
Agaev:2020zad .

The correlator $\Pi^{\mathrm{OPE}}(p)$ has also a simple structure
$\sim\mathrm{I}$ and is characterized by an amplitude
$\Pi^{\mathrm{OPE}}(p^{2})$. To find a preliminary sum rule, we equate the
amplitudes $\Pi^{\mathrm{Phys}}(p^{2})$ and $\Pi^{\mathrm{OPE}}(p^{2})$. This
equality contains contributions coming from both the ground-state particle and
higher resonances. The latter can be suppressed by applying the Borel
transformation to both sides of the sum rule equality. Afterward, using the
quark-hadron duality assumption, we subtract the suppressed terms from the
obtained expression. These manipulations lead to dependence of the sum rule
equality on the Borel and continuum subtraction (threshold) parameters $M^{2}$
and $s_{0}$. Obtained by this way, the expression and its derivative over
$d/d(-1/M^{2})$ allow us to find the sum rules for the mass $m$ and coupling
$f$ of the molecule $M$, which 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)

The function $\Pi(M^{2},s_{0})$ in Eqs. (11) and (12) is the invariant
amplitude $\Pi^{\mathrm{OPE}}(p^{2})$ after the Borel transformation and
continuum subtraction, and
$\Pi^{\prime}(M^{2},s_{0})=d\Pi(M^{2},s_{0})/d(-1/M^{2})$.

The Borel transform of the amplitude $\Pi^{\mathrm{Phys}}(p^{2})$ is given by
the formula

$\mathcal{B}\Pi^{\mathrm{Phys}}(p^{2})=fme^{-m^{2}/M^{2}}.$ (13)

For the correlator $\Pi(M^{2},s_{0})$, we find a more complicated expression

$\Pi(M^{2},s_{0})=\int_{\mathcal{M}^{2}}^{s_{0}}ds\rho^{\mathrm{OPE}}(s)e^{-s/M^{2}}+\Pi(M^{2}),$
(14)

where $\mathcal{M}=m_{c}+m_{s}$ is the mass of constituent quarks in the
molecule $M$. It is worth noting, that we neglect masses of $u$ and $d$
quarks, but take into account terms $\sim m_{s}$. At the same time, we do not
include into analysis contributions $\sim m_{s}^{2}$ and set $m_{s}^{2}=0$.
The spectral density $\rho^{\mathrm{OPE}}(s)$ is computed as an imaginary part
of the amplitude $\Pi^{\mathrm{OPE}}(p^{2})$. Borel transforms some of terms
obtained directly from $\Pi^{\mathrm{OPE}}(p)$ are denoted in Eq. (14) by
$\Pi(M^{2}).$ In this paper, calculations are carried out by taking into
account vacuum condensates up to dimension $10$. Analytical expressions of
$\rho^{\mathrm{OPE}}(s)$ and $\Pi(M^{2})$ are lengthy, therefore we do not
write down them here explicitly.

For numerical computations of $m$ and $f$, one should specify different vacuum
condensates, which enter to the sum rules in Eqs. (11) and (12). These
condensates are universal parameters, which were extracted from analysis of
numerous processes: Their numerical values 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{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\langle
g_{s}^{3}G^{3}\rangle=(0.57\pm 0.29)~{}\mathrm{GeV}^{6},$ $\displaystyle
m_{s}=93_{-5}^{+11}~{}\mathrm{MeV},\ m_{c}=1.27\pm 0.02~{}\mathrm{GeV}.$ (15)

We have also included into Eq. (15) the masses of $c$ and $s$ quarks.

Figure 1: Dependence of the pole contribution on the Borel parameter $M^{2}$
at fixed $s_{0}$. The limit $\mathrm{PC}=0.5$ is shown by the horizontal black
line. The red diamond notes the point, where the mass $m$ of the molecule
$M=D^{\ast+}K^{\ast+}$ has effectively been calculated.

The sum rules in Eqs. (11) and (12) depend also on the Borel and continuum
threshold parameters $M^{2}$ and $s_{0}$. The choice of working windows for
$M^{2}$ and $s_{0}$ has to satisfy standard constraints imposed on the pole
contribution ($\mathrm{PC}$) and convergence of the operator product
expansion. To quantify these constraints, it is appropriate to use expressions

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

and

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

where $\Pi^{\mathrm{DimN}}(M^{2},s_{0})=$
$\Pi^{\mathrm{Dim(8+9+10)}}(M^{2},s_{0})$.

The $\mathrm{PC}$ is used to fix maximum of the Borel parameter
$M_{\mathrm{max}}^{2}$, whereas its minimal value $M_{\mathrm{min}}^{2}$ is
limited by the convergence of $\mathrm{OPE}$. In sum rule analyses of ordinary
hadrons $\mathrm{PC}\geq 0.5$ is a standard requirement. When studying
multiquark particles this constraint reduces the region of the allowed
$M^{2}$. To be convinced in convergence of the operator product expansion, we
demand fulfillment of $R(M_{\mathrm{min}}^{2})\leq 0.05$.

Our calculations prove that the regions for the parameters $M^{2}$ and
$s_{0}$,

$M^{2}\in[2,2.7]~{}\mathrm{GeV}^{2},\ s_{0}\in[10.7,11.7]~{}\mathrm{GeV}^{2},$
(18)

meet all required restrictions. Thus, at $M^{2}=2.7~{}\mathrm{GeV}^{2}$ the
pole contribution on average in $s_{0}$ is $0.53$, whereas at
$M^{2}=2~{}\mathrm{GeV}^{2}$ it equals to $0.72$. In Fig. 1 the pole
contribution is plotted as a function of $M^{2}$ at various fixed $s_{0}$. It
is seen that by excluding the small region $M^{2}\geq 2.65~{}\mathrm{GeV}^{2}$
at $s_{0}=10.7~{}\mathrm{GeV}^{2}$ the pole contribution exceeds $0.5$. On
average in $s_{0}$, the constraint $\mathrm{PC}\geq 0.5$ is satisfied in the
all working window for the Borel parameter. At the minimum point
$M^{2}=2~{}\mathrm{GeV}^{2}$, we find $R(2~{}\mathrm{GeV}^{2})\approx 0.015$
and the sum of dimension-$8,9,10$ contributions is less than $1.5\%$ of the
full result.

Dominance of the perturbative contribution to $\Pi(M^{2},s_{0})$, and
convergence of the operator product expansion are another important problems
in the sum rule studies. In Fig. 2, we compare the perturbative and
nonperturbative components of the correlation function. One sees that the
perturbative contribution to $\Pi(M^{2},s_{0})$ prevails over nonperturbaative
one, and forms more than $53\%$ of $\Pi(M^{2},s_{0})$ already at
$M^{2}=2~{}\mathrm{GeV}^{2}$ growing gradually in the considered range of
$M^{2}$. From this figure it is also clear that convergence of $\mathrm{OPE}$
is satisfied: Contributions of the nonperturbative terms reduce by increasing
the dimensions of the corresponding operators. There is some disordering in
these contributions connected with smallness of gluon condensates. The
dimension-$3,\ 6,\ 9$ and $10$ terms are positive. The $\mathrm{Dim3}$ and
$\mathrm{Dim6}$ contributions numerically exceed contributions of other
operators, whereas $\mathrm{Dim9}$ and $\mathrm{Dim10}$ terms are very small
and not shown in the plot.

Predictions for $m$ and $f$ are obtained by taking the mean values of these
parameters calculated at different choices of $M^{2}$ and $s_{0}$:

$\displaystyle m$ $\displaystyle=$ $\displaystyle(2924~{}\pm
107)~{}\mathrm{MeV},$ $\displaystyle f$ $\displaystyle=$ $\displaystyle(4.3\pm
0.7)\times 10^{-3}~{}\mathrm{GeV}^{4}.$ (19)

The $m$ and $f$ from Eq. (19) effectively correspond to the sum rules’ results
at $M^{2}=2.3~{}\mathrm{GeV}^{2}$ and $s_{0}=11.2~{}\mathrm{GeV}^{2}$ noted in
Fig. 1 by the red diamond. This point is approximately at the middle of the
intervals shown in Eq. (18), where the pole contribution is
$\mathrm{PC}\approx 0.62$. The circumstances discussed above ensure the
ground-state nature of $M$ and reliability of the obtained results.

The dependence of the mass $m$ on the parameters $M^{2}$ and $s_{0}$ is drawn
in Fig. 3. In general, a physical quantity should not depend on the auxiliary
parameter of computations $M^{2}$. Nevertheless, such residual dependence of
$m$ on the Borel parameter, as well as on $s_{0}$ exists and generates
theoretical ambiguities of the extracted predictions in Eq. (19). It is worth
noting that these ambiguities are smaller for $m$ than for $f$. Indeed, for
the mass they are equal only to $\pm 4\%$ of the central value, whereas in the
case of $f$, they amount to $\pm 16\%$. Such difference is connected by the
analytical forms of the sum rules for these quantities: The mass $m$ is given
by the ratio of the correlation functions which smooths relevant effect,
whereas $f$ depends on $\Pi(M^{2},s_{0})$ itself.

Our result for the mass of the molecule $D^{\ast+}K^{\ast+}$ agrees very well
with the LHCb datum. This is necessary, but not enough to make credible
conclusions about the nature of the resonance $T_{cs0}^{a++}$: For more
reliable statements, one needs to estimate also the full width of the hadronic
molecule $D^{\ast+}K^{\ast+}$ suggested in this paper to model
$T_{cs0}^{a++}$.

Figure 2: Left: Different contributions to $\Pi(M^{2},s_{0})$ normalized to
$1$ as functions of the Borel parameter $M^{2}$, Right: Normalized
contributions of various operators to $\Pi(M^{2},s_{0})$ as functions of
$M^{2}$. All curves in this figure have been calculated at
$s_{0}=11.2~{}\mathrm{GeV}^{2}$.

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

##  III Width of the molecule $M$

The quark content, spin-parity and mass of the molecule $M=$
$D^{\ast+}K^{\ast+}$ allow us to classify its decay channels. Because
$T_{cs0}^{a++}$ was seen in the $D_{s}^{+}\pi^{+}$ invariant mass
distribution, the process $M\rightarrow D_{s}^{+}\pi^{+}$ should be the
dominant decay channel of the molecule $M$. The processes $M\rightarrow
D_{s}^{\ast+}\rho^{+}$ and $M\rightarrow D^{\ast+}K^{\ast+}$ are also among
the kinematically allowed decay modes of the $M$. We are going to evaluate the
full width of the molecule $M$ using these decay channels. It is noting that,
$\rho^{+}$ and $K^{\ast+}$ mesons decay almost exclusively to $\pi^{+}\pi^{0}$
and $(K\pi)^{+}$ pairs (see, Ref. PDG:2022 ), therefore widths of the last two
processes can be considered also as widths of the modes $M\rightarrow
D_{s}^{\ast+}\pi^{+}\pi^{0}$ and $M\rightarrow D^{\ast+}(K\pi)^{+}$,
respectively.

Partial widths of aforementioned processes are determined by the strong
couplings $g_{i}$ at vertices $MD_{s}^{+}\pi^{+}$, and etc. One of the
effective ways to evaluate them is the QCD light-cone sum rule method
Balitsky:1989ry . In this approach, the QCD side of the sum rule, instead of
the local vacuum condensates, is expressed in terms of one of the final
mesons’ distribution amplitudes (DAs). In general, the DAs of a hadron are
nonlocal matrix elements of various operators with different twists sandwiched
between the hadron and vacuum states. They are specific for each particle and
modeled by taking into account the available experimental data. The LCSR
method is suitable for analysis of not only the conventional hadrons, but also
the multiquark systems such as tetraquarks, pentaquarks, etc. In the case of
tetraquark-meson-meson vertices, however, these DAs reduce to the local matrix
elements of the meson. For treatment of such vertices, one needs to apply some
additional mathematical tools. The LCSR method was adapted for investigation
of the tetraquark-meson-meson vertices in Ref. Agaev:2016dev  and applied to
study numerous decays of four-quark exotic mesons Agaev:2020zad .

###  III.1 Decay $M\rightarrow D_{s}^{+}\pi^{+}$

Here, we consider, in a detailed manner, the decay of the molecule $M$ to pair
of pseudoscalar mesons $D_{s}^{+}\pi^{+}$. Partial width of this process
depends on the spectroscopic parameters of the initial and final-state
particles. The spectroscopic parameters of $M$ have been evaluated in the
previous section. The masses and decay constants of the mesons $D_{s}^{+}$ and
$\pi^{+}$ are available from other sources. The only unknown quantity required
to calculate the $M\rightarrow D_{s}^{+}\pi^{+}$ decays’s partial width is the
strong coupling $g_{1}$ of the particles at the vertex $MD_{s}^{+}\pi^{+}$.

The coupling $g_{1}$ is defined in terms of the on-mass-shell matrix element,

$\langle\pi(q)D_{s}\left(p\right)|M(p^{\prime})\rangle=g_{1}p\cdot
p^{\prime},$ (20)

where the mesons $\pi^{+}$ and $D_{s}^{+}$ are denoted as $\pi$ and $D_{s}$,
respectively. In Eq. (20) $p^{\prime}$, $p$ and $q$ are four-momenta of $M$,
and mesons $D_{s}$ and $\pi$.

In the framework of the LCSR method the sum rule for the coupling $g_{1}$ can
be obtained from the correlation function,

$\Pi(p,q)=i\int
d^{4}xe^{ipx}\langle\pi(q)|\mathcal{T}\\{J^{D_{s}}(x)J^{{\dagger}}(0)\\}|0\rangle,$
(21)

where $J(x)$ is the current for the molecule $M$ given by Eq. (2). The
interpolating current $J^{D_{s}}(x)$ for the meson $D_{s}^{+}$ is defined by
the formula

$J^{D_{s}}(x)=\overline{s}_{j}(x)i\gamma_{5}c_{j}(x),$ (22)

with $j$ being the color index.

At this stage of our analysis, we express the correlation function $\Pi(p,q)$
using the physical parameters of the particles involved into the decay
process. To this end, we write $\Pi^{\mathrm{Phys}}(p,q)$ in the factorized
form Belyaev:1994zk ; Agaev:2022iha ,

$\displaystyle\Pi^{\mathrm{Phys}}(p,q)$ $\displaystyle=$
$\displaystyle\frac{\langle\pi(q)D_{s}\left(p\right)|M(p^{\prime})\rangle\langle
M(p^{\prime})|J^{\dagger}|0\rangle}{\left(p^{\prime 2}-m^{2}\right)}$ (23)
$\displaystyle\times\frac{\ \langle
0|J^{D_{s}}|D_{s}(p)\rangle}{\left(p^{2}-m_{D_{s}}^{2}\right)},$

where $m_{D_{s}}$ is the mass of the meson $D_{s}^{+}$. As is seen, the
function $\Pi^{\mathrm{Phys}}(p,q)$ contains the matrix elements of the vertex
$MD_{s}^{+}\pi^{+}$, the molecule $M$ and meson $D_{s}^{+}$. The matrix
element of $M$ is known from Eq. (8), whereas for the meson $D_{s}^{+}$ we use

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

with $f_{D_{s}}$ being the decay constant of $D_{s}^{+}$.

After simple manipulations, we get

$\displaystyle\Pi^{\mathrm{Phys}}(p,q)=g_{1}\frac{fmf_{D_{s}}m_{D_{s}}^{2}}{(m_{c}+m_{s})\left(p^{2}-m_{D_{s}}^{2}\right)}$
$\displaystyle\times\frac{1}{\left(p^{\prime 2}-m^{2}\right)}p\cdot
p^{\prime}+\cdots.$ (25)

The term in Eq. (25) corresponds to the contribution of ground-state particles
in $M$ and $D_{s}^{+}$ channels. Effects of higher resonances and continuum
states in these channels are denoted by ellipses. The function
$\Pi^{\mathrm{Phys}}(p,q)$ form the physical side of the sum rule for the
coupling $g_{1}$. It has the Lorentz structure proportional to the unit matrix
$\mathrm{I}$, therefore the expression in r.h.s. of Eq. (25) is the invariant
amplitude $\Pi^{\mathrm{Phys}}(p^{2},p^{\prime 2})$, which depends on two
variables $p^{2}$ and $p^{\prime 2}$.

The correlation function $\Pi(p,q)$ calculated in terms of quark propagators
and matrix elements of the pion constitutes the QCD side of sum rules and is
equal to

$\displaystyle\Pi^{\mathrm{OPE}}(p,q)=-\int
d^{4}xe^{ipx}\left[\gamma_{\mu}S_{s}^{bj}(-x){}\gamma_{5}S_{c}^{ja}(x)\gamma^{\mu}\right]_{\alpha\beta}$
$\displaystyle\times\langle\pi(q)|\overline{u}_{\alpha}^{b}(0)d_{\beta}^{a}(0)|0\rangle,$
(26)

where $\alpha$ and $\beta$ are the spinor indices.

The correlator $\Pi^{\mathrm{OPE}}(p,q)$, apart from propagators of the $c$
and $s$ quarks, contains also the local matrix elements
$\langle\pi|\overline{u}_{\alpha}^{b}d_{\beta}^{a}|0\rangle$ of the pion. In
the standard LCSR method, while studying vertices of conventional mesons, the
correlator depends on the nonlocal matrix elements of the meson (for instance,
$\pi$ meson), which after some transformations can be expressed in terms of
its different DAs. In the case under discussion, the $\Pi^{\mathrm{OPE}}(p,q)$
contains the pion’s local matrix elements appearance of which has a simple
explanation. Indeed, because the hadronic molecule $M$ is built of four
valence quarks located at space-time position $x=0$, contractions of two quark
operators from the currents $J^{D_{s}}(x)$ and $J(x)$ leave free the two quark
fields from $M$ at the same position $x=0$. This feature of the correlation
function $\Pi^{\mathrm{OPE}}(p,q)$ connected with differences in the quark
contents of tetraquarks and ordinary mesons is unavoidable effect for all
tetraquark-meson-meson vertices.

It turns out that the $\Pi^{\mathrm{OPE}}(p,q)$-type correlators emerge in the
limit $q\rightarrow 0$ in LCSR calculations Belyaev:1994zk , which is known as
the soft-meson approximation. In this approximation, $p=p^{\prime}$ and
invariant amplitudes $\Pi^{\mathrm{Phys}}(p^{2})$ and
$\Pi^{\mathrm{OPE}}(p^{2})$ depend only on one variable $p^{2}$. It is worth
emphasizing that the limit $q\rightarrow 0$ is applied to hard parts of the
invariant amplitudes, whereas in their soft parts (i.e., in matrix elements)
terms $\sim q^{2}=m^{2}$ are taken into account. In other words, soft-meson
approximation should not be considered as a massless limit of the correlation
functions.

Technical difficulties generated by this limit in the physical side of sum
rules can be cured by means of technical tools elaborated in Refs.
Ioffe:1983ju ; Belyaev:1994zk . It is important that the sum rules for the
strong couplings obtained using the full LCSR method, and soft-meson
approximation lead to numerically close results Belyaev:1994zk . To clarify
this last point, we note that in the full version of the LCSR method a sum
rule for the strong coupling at a vertex of three conventional mesons depends
on numerous two- and three-particle quark-gluon DAs of a final meson. In the
limit $q\rightarrow 0$ in this expression survive only a few leading terms.
Because their contributions are numerically decisive, for the strong coupling
the full version and soft-meson limit of the light-cone sum rules give close
predictions. Thus, in Ref. Belyaev:1994zk  the couplings $g_{D^{\ast}D\pi}$
and $g_{B^{\ast}B\pi}$ at the vertices $D^{\ast}D\pi$ and $B^{\ast}B\pi$ were
calculated using both of these methods. In the full LCSR approach these
couplings are equal to

$g_{D^{\ast}D\pi}^{\mathrm{full}}=12.5\pm 1.0,\
g_{B^{\ast}B\pi}^{\mathrm{full}}=29\pm 3,$ (27)

whereas in the soft-meson approximation, the authors found

$g_{D^{\ast}D\pi}^{\mathrm{soft}}=11\pm 2,\
g_{B^{\ast}B\pi}^{\mathrm{soft}}=28\pm 6.$ (28)

As is seen, values of the couplings are very close to each other, though
uncertainties in the soft-limit are larger than in the full version
Belyaev:1994zk .

These arguments do not imply a necessity of the soft-meson approximation to
study all of the vertices containing four-quark states. Two tetraquark-meson
vertices, for example, can be readily investigated in the context of the
conventional LCSR method Agaev:2016srl .

The local matrix elements,
$\langle\pi|\overline{u}_{\alpha}^{b}d_{\beta}^{a}|0\rangle$, carry color and
spinor indices and are uneasy objects for further operations. We can rewrite
them in convenient forms by expanding $\overline{u}d$ over the full set of
Dirac matrices $\Gamma^{J}$,

$\Gamma^{J}=\mathbf{1},\ \gamma_{5},\ \gamma_{\mu},\ i\gamma_{5}\gamma_{\mu},\
\sigma_{\mu\nu}/\sqrt{2},$ (29)

and projecting onto the colorless states

$\overline{u}_{\alpha}^{b}(0)d_{\beta}^{a}(0)\rightarrow\frac{1}{12}\delta^{ba}\Gamma_{\beta\alpha}^{J}\left[\overline{u}(0)\Gamma^{J}d(0)\right].$
(30)

The operators $\overline{u}\Gamma^{J}d$, sandwiched between the $\pi$ meson
and vacuum, generate local matrix elements of the $\pi$ meson, which can be
implemented into $\Pi^{\mathrm{OPE}}(p,q)$.

The expression obtained for the $\Pi^{\mathrm{OPE}}(p^{2})$ in the soft limit
is considerably more simple than the invariant amplitude in the full version
of LCSR method. But, at the same time, soft-meson approximation produces
problems in the physical side of the sum rule. Below, we will come back to
finish calculation of $\Pi^{\mathrm{OPE}}(p^{2})$, but now it is time to fix
sources of complications in $\Pi^{\mathrm{Phys}}(p^{2},p^{\prime 2})$. To this
end, we rewrite this amplitude in the soft limit

$\displaystyle\Pi^{\mathrm{Phys}}(p^{2})=g_{1}\frac{fmf_{D_{s}}m_{D_{s}}^{2}}{2(m_{c}+m_{s})\left(p^{2}-\widetilde{m}^{2}\right)^{2}}$
$\displaystyle\times(2\widetilde{m}^{2}-m_{\pi}^{2})+\cdots,$ (31)

where $\widetilde{m}^{2}=(m^{2}+m_{D_{s}}^{2})/2$ and $m_{\pi}$ is the mass of
the pion. It is seen that the amplitude $\Pi^{\mathrm{Phys}}(p^{2})$, in the
soft approximation, instead of two poles at different points has one double
pole at $p^{2}=\widetilde{m}^{2}$.

The Borel transformation of $\Pi^{\mathrm{Phys}}(p^{2})$ is given by the
expression,

$\displaystyle\Pi^{\mathrm{Phys}}(M^{2})=\left[g_{1}\frac{fmf_{D_{s}}m_{D_{s}}^{2}}{2(m_{c}+m_{s})}(2\widetilde{m}^{2}-m_{\pi}^{2})\right.$
$\displaystyle\left.+AM^{2}\right]\frac{e^{-\widetilde{m}^{2}/M^{2}}}{M^{2}}+\cdots.$
(32)

Apart from the ground-state contribution, in the soft limit, the amplitude
$\Pi^{\mathrm{Phys}}(M^{2})$ contains additional unsuppressed terms $\sim A$
Belyaev:1994zk . These terms correspond to the transitions from the excited
states in $M=D^{\ast+}K^{\ast+}$ channel with $m^{\ast}>m$, and are not
suppressed relative to the ground-state contribution even after the Borel
transformation. These circumstances make problematic extraction of the $g_{1}$
from Eq. (32). The contributions $\sim A$ can be removed from the physical
side of the sum rule by means of the operator
$\mathcal{P}(M^{2},\widetilde{m}^{2})$ Ioffe:1983ju ; Belyaev:1994zk ,

$\mathcal{P}(M^{2},\widetilde{m}^{2})=\left(1-M^{2}\frac{d}{dM^{2}}\right)M^{2}e^{\widetilde{m}^{2}/M^{2}},$
(33)

which should be applied to both sides of the sum rule equality.

After this operation the remaining suppressed terms in
$\Pi^{\mathrm{Phys}}(M^{2})$, denoted in Eq. (32) by ellipses, can be
subtracted by a standard way. As a result, we find the sum rule for the strong
coupling $g_{1}$, which reads

$g_{1}=\frac{2(m_{c}+m_{s})}{fmf_{D_{s}}m_{D_{s}}^{2}(2\widetilde{m}^{2}-m_{\pi}^{2})}\mathcal{P}(M^{2},\widetilde{m}^{2})\Pi^{\mathrm{OPE}}(M^{2},s_{0}),$
(34)

where $\Pi^{\mathrm{OPE}}(M^{2},s_{0})$ is the Borel transformed and
subtracted invariant amplitude $\Pi^{\mathrm{OPE}}(p^{2})$.

Recipes to compute the correlation function $\Pi^{\mathrm{OPE}}(p,q)$ in the
soft approximation were explained in Ref. Agaev:2016dev , therefore we give
only principal points of these calculations. Thus, having substituted the
expansion (30) into Eq. (26), we perform summations over color indices and fix
the local matrix elements of the pion that contribute to the
$\Pi^{\mathrm{OPE}}(p,q)$ in the soft-meson limit. It turns out that the
contributions to $\Pi^{\mathrm{OPE}}(p,q=0)$ come from the pion’s two-particle
twist-3 matrix element,

$\langle 0|\overline{d}i\gamma_{5}u|\pi\rangle=f_{\pi}\mu_{\pi},$ (35)

where

$\mu_{\pi}=\frac{m_{\pi}^{2}}{m_{u}+m_{d}}=-\frac{2\langle\overline{q}q\rangle}{f_{\pi}^{2}}.$
(36)

The second equality in Eq. (36) arises from the partial conservation of axial-
vector current (PCAC).

The amplitude $\Pi^{\mathrm{OPE}}(M^{2},s_{0})$ is given by the formula

$\displaystyle\Pi^{\mathrm{OPE}}(M^{2},s_{0})=\frac{f_{\pi}\mu_{\pi}}{8\pi^{2}}\int_{\mathcal{M}^{2}}^{s_{0}}\frac{ds(m_{c}^{2}-s)}{s}$
$\displaystyle\times(m_{c}^{2}-2m_{c}m_{s}-s)e^{-s/M^{2}}+f_{\pi}\mu_{\pi}\Pi_{\mathrm{NP}}(M^{2}).$
(37)

The first term in Eq. (37) is the perturbative component of the
$\Pi^{\mathrm{OPE}}(M^{2},s_{0})$. The nonperturbative term
$\Pi_{\mathrm{NP}}(M^{2})$ is computed with dimension-$9$ accuracy and has the
following form:

$\displaystyle\Pi_{\mathrm{NP}}(M^{2})=\frac{\langle\overline{s}s\rangle}{6M^{2}}\left[M^{2}(2m_{c}-m_{s})-m_{c}^{2}m_{s}\right]e^{-m_{c}^{2}/M^{2}}$
$\displaystyle+\langle\frac{\alpha_{s}G^{2}}{\pi}\rangle\frac{m_{c}}{72M^{4}}\int_{0}^{1}\frac{dx}{(1-x)^{3}}\left[m_{s}(1-x)(2M^{2}+m_{c}^{2})\right.$
$\displaystyle\left.-m_{c}^{3}x\right]e^{-m_{c}^{2}/[M^{2}(1-x)]}-\frac{\langle\overline{s}g\sigma
Gs\rangle m_{c}^{3}}{6M^{6}}(M^{2}-m_{c}m_{s})$ $\displaystyle\times
e^{-m_{c}^{2}/M^{2}}+\langle\frac{\alpha_{s}G^{2}}{\pi}\rangle\langle\overline{s}s\rangle\frac{\pi^{2}m_{c}}{9M^{8}}\left[M^{2}m_{c}(m_{c}+m_{s})\right.$
$\displaystyle\left.-2m_{c}^{3}m_{s}-M^{4}\right]e^{-m_{c}^{2}/M^{2}}+\langle\frac{\alpha_{s}G^{2}}{\pi}\rangle\langle\overline{s}g\sigma
Gs\rangle\frac{\pi^{2}m_{c}}{18M^{12}}$
$\displaystyle\times\left[3M^{6}-20m_{c}^{5}m_{s}+10m_{c}^{3}M^{2}(m_{c}+2m_{s})\right.$
$\displaystyle\left.-4M^{4}m_{c}(3m_{c}+m_{s})\right]e^{-m_{c}^{2}/M^{2}}.$
(38)

Parameters | Values (in $\mathrm{MeV}$ units)
---|---
$m_{D_{s}}$ | $1969.0\pm 1.4$
$m_{\pi}$ | $139.57039\pm 0.00017$
$m_{D_{s}^{\ast}}$ | $2112.2\pm 0.4$
$m_{\rho}$ | $775.11\pm 0.34$
$m_{D^{\ast}}$ | $2010.26\pm 0.05$
$m_{K^{\ast}}$ | $891.67\pm 0.26$
$f_{D_{s}}$ | $249.9\pm 0.5$
$f_{\pi}$ | $130.2\pm 0.8$
$f_{D_{s}^{\ast}}$ | $268.8\pm 6.5$
$f_{\rho}$ | $216\pm 3$
$f_{D^{\ast}}$ | $223.5\pm 0.5$
$f_{K^{\ast}}$ | $204\pm 0.3$
Table 1: Masses and decay constants of the mesons $D_{s}^{(\ast)}$,
$D^{\ast}$, $K^{\ast}$, $\pi$ and $\rho$, which have been used in numerical
computations.

Besides the vacuum condensates, the sum rule in Eq. (34) contains also the
masses and decay constants of the final-state mesons $D_{s}^{+}$ and
$\pi^{+}$. Numerical values of these parameters, as well as parameters of
other mesons, which will be necessary later, are presented in Table 1. For the
decay constants of the mesons $D_{s}^{\ast+}$ and $D^{\ast+}$, we employ
predictions of the QCD lattice method Lubicz:2016bbi . For all other
parameters, we use information of the Particle Data Group, mainly its last
edition Ref. PDG:2022 . All these input parameters were either measured
experimentally or extracted by means of alternative theoretical approaches.

Numerical analysis demonstrates that the working regions shown in Eq. (18)
used in calculations of the $M$ molecule’s mass meet the required restrictions
on the Borel and continuum subtraction parameters $M^{2}$ and $s_{0}$ imposed
in the case of the decay process. Therefore, in computations of
$\Pi^{\mathrm{OPE}}(M^{2},s_{0})$ the parameters $M^{2}$ and $s_{0}$ have been
varied within the limits (18).

For $g_{1}$, the numerical calculations yield

$g_{1}=(7.1\pm 1.1)\times 10^{-1}~{}\mathrm{GeV}^{-1}.$ (39)

The partial width of the decay $M\rightarrow D_{s}^{+}\pi^{+}$ is determined
by the expression

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

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

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

Then it is not difficult to find that

$\Gamma_{1}\left[M\rightarrow D_{s}^{+}\pi^{+}\right]=(71\pm
23)~{}\mathrm{MeV},$ (42)

which is large enough to confirm dominant nature of this channel.

###  III.2 Processes $M\rightarrow D_{s}^{\ast+}\rho^{+}$ and $M\rightarrow
D^{\ast+}K^{\ast+}$

These two processes differ from previous decay by a vector nature of produced
mesons. This fact modifies matrix elements of the vertices and formulas for
decay widths of the processes. We concentrate on analysis of the decay
$M\rightarrow D_{s}^{\ast+}\rho^{+}$, and write down only the final
predictions for $M\rightarrow D^{\ast+}K^{\ast+}$.

The correlation function, which should be studied in the case of the decay
$M\rightarrow D_{s}^{\ast+}\rho^{+}$, has the following form

$\Pi_{\mu}(p,q)=i\int
d^{4}xe^{ipx}\langle\rho(q)|\mathcal{T}\\{J_{\mu}^{D_{s}^{\ast}}(x)J^{{\dagger}}(0)\\}|0\rangle,$
(43)

where $J_{\mu}^{D_{s}^{\ast}}(x)$ is the interpolating current for the meson
$D_{s}^{\ast+}$

$J_{\mu}^{D_{s}^{\ast}}(x)=\overline{s}_{i}(x)\gamma_{\mu}c_{i}(x).$ (44)

The correlation function $\Pi_{\mu}(p,q)$ written down in terms of the matrix
elements of the particles and the vertex $MD_{s}^{\ast+}\rho^{+}$ is given by
the expression

$\displaystyle\Pi_{\mu}^{\mathrm{Phys}}(p,q)=\frac{\langle\rho(q,\epsilon)D_{s}^{\ast+}(p,\varepsilon)|M(p^{\prime})\rangle\langle
M(p^{\prime})|J^{\dagger}|0\rangle}{(p^{\prime 2}-m^{2})}$
$\displaystyle\times\frac{\langle
0|J_{\mu}^{D_{s}^{\ast}}|D_{s}^{\ast+}(p,\varepsilon)\rangle}{(p^{2}-m_{D_{s}^{\ast}}^{2})}+\cdots,$
(45)

where $m_{D_{s}^{\ast}}$ and $\varepsilon_{\mu}$ are the mass and polarization
vector of the meson $D_{s}^{\ast+}$, and $\epsilon_{\nu}$ is polarization
vector of the $\rho$ meson. The expression in Eq. (45) is the contribution of
the ground-state particles to the physical side of the sum rule, whereas
ellipses stand for the contributions of higher resonances and continuum
states.

We introduce the matrix element of the meson $D_{s}^{\ast+}$ by the formula

$\langle
0|J_{\mu}^{D_{s}^{\ast}}|D_{s}^{\ast+}(p,\varepsilon)\rangle=m_{D_{s}^{\ast}}f_{D_{s}^{\ast}}\varepsilon_{\mu}.$
(46)

We also model the mass-shell matrix element of the vertex
$MD_{s}^{\ast+}\rho^{+}$ in the following way:

$\displaystyle\langle\rho(q,\epsilon)D_{s}^{\ast+}(p,\varepsilon)|M(p^{\prime})\rangle=g_{2}\left[\left(q\cdot
p\right)\left(\epsilon\cdot\varepsilon^{\ast}\right)\right.$
$\displaystyle\left.-\left(p\cdot\epsilon\right)\left(q\cdot\varepsilon^{\ast}\right)\right],$
(47)

Then, it is not difficult to calculate the function
$\Pi_{\mu}^{\mathrm{Phys}}(p,q)$,

$\displaystyle\Pi_{\mu}^{\mathrm{Phys}}(p,q)=\frac{fmm_{D_{s}^{\ast}}f_{D_{s}^{\ast}}}{(p^{\prime
2}-m^{2})(p^{2}-m_{D_{s}^{\ast}}^{2})}$
$\displaystyle\times\left(\frac{m_{\rho}^{2}+m_{D_{s}^{\ast}}^{2}-m^{2}}{2}\epsilon_{\mu}+p\cdot\epsilon
q_{\mu}\right)+\cdots,$ (48)

As is seen, $\Pi_{\mu}^{\mathrm{Phys}}(p,q)$ contains two structures
proportional to $\epsilon_{\mu}$ and $q_{\mu}$. We are going to employ the
structure $\sim\epsilon_{\mu}$, and corresponding invariant amplitude which in
the soft limit has the form

$\widehat{\Pi}^{\mathrm{Phys}}(p^{2})=\frac{fmm_{D_{s}^{\ast}}f_{D_{s}^{\ast}}}{(p^{2}-\widehat{m}^{2})^{2}}\frac{m_{\rho}^{2}+m_{D_{s}^{\ast}}^{2}-m^{2}}{2},$
(49)

with $\widehat{m}^{2}$ being equal to $(m^{2}+m_{D_{s}^{\ast}}^{2})/2$.

The QCD side of the sum rule is determined by the correlator

$\displaystyle\Pi_{\mu}^{\mathrm{OPE}}(p,q)=i\int
d^{4}xe^{ipx}\left[\gamma_{\nu}S_{s}^{bi}(-x){}\gamma_{\mu}S_{c}^{ia}(x)\gamma^{\nu}\right]_{\alpha\beta}$
$\displaystyle\times\langle\rho(q)|\overline{u}_{\alpha}^{b}(0)d_{\beta}^{a}(0)|0\rangle.$
(50)

This function has the same Lorentz structures as
$\Pi_{\mu}^{\mathrm{Phys}}(p,q)$. In the soft-meson approximation, it receives
contribution from the matrix element

$\langle
0|\overline{d}\gamma_{\nu}u|\rho\rangle=f_{\rho}m_{\rho}\epsilon_{\nu},$ (51)

where $f_{\rho}$ is the decay constant of the $\rho$ meson.

Having fixed an amplitude which is proportional to $\epsilon_{\mu}$ and
labeled it by $\widehat{\Pi}^{\mathrm{OPE}}(p^{2})$ it is not difficult to
write the sum rule for the strong coupling $g_{2}$

$g_{2}=\frac{2\mathcal{P}(M^{2},\widehat{m}^{2})\widehat{\Pi}^{\mathrm{OPE}}(M^{2},s_{0})}{fmm_{D_{s}^{\ast}}f_{D_{s}^{\ast}}(m_{\rho}^{2}+m_{D_{s}^{\ast}}^{2}-m^{2})}.$
(52)

Here, $\widehat{\Pi}^{\mathrm{OPE}}(M^{2},s_{0})$ is the amplitude
$\widehat{\Pi}^{\mathrm{OPE}}(p^{2})$ after the Borel transformation and
continuum subtraction procedures. It does not differ considerably from the
$\Pi^{\mathrm{OPE}}(M^{2},s_{0})$ and has the following form

$\displaystyle\widehat{\Pi}^{\mathrm{OPE}}(M^{2},s_{0})=\frac{f_{\rho}m_{\rho}}{48\pi^{2}}\int_{\mathcal{M}^{2}}^{s_{0}}\frac{ds(m_{c}^{2}-s)}{s^{2}}$
$\displaystyle\times(m_{c}^{4}+m_{c}^{2}m_{s}-6m_{c}m_{s}s-2s^{2})e^{-s/M^{2}}$
$\displaystyle+f_{\rho}m_{\rho}\widehat{\Pi}_{\mathrm{NP}}(M^{2}).$ (53)

The nonperturbative term $\widehat{\Pi}_{\mathrm{NP}}(M^{2})$ is given by the
expression

$\displaystyle\widehat{\Pi}_{\mathrm{NP}}(M^{2})=\frac{\langle\overline{s}s\rangle
m_{c}}{12M^{2}}\left(2M^{2}-m_{c}m_{s}\right)e^{-m_{c}^{2}/M^{2}}$
$\displaystyle+\langle\frac{\alpha_{s}G^{2}}{\pi}\rangle\frac{m_{c}}{144M^{4}}\int_{0}^{1}\frac{dx}{(1-x)^{3}}\left[m_{c}^{3}x+(1-x)\left(m_{c}^{2}m_{s}\right.\right.$
$\displaystyle\left.\left.-2m_{s}M^{2}-m_{c}M^{2}\right)\right]e^{-m_{c}^{2}/[M^{2}(1-x)]}+\frac{\langle\overline{s}g\sigma
Gs\rangle}{72M^{6}}$
$\displaystyle\times(6m_{c}^{4}m_{3}-6m_{c}^{3}M^{2}-2m_{c}^{2}m_{s}M^{2}-m_{s}M^{4})e^{-m_{c}^{2}/M^{2}}$
$\displaystyle+\langle\frac{\alpha_{s}G^{2}}{\pi}\rangle\langle\overline{s}s\rangle\frac{\pi^{2}m_{c}}{36M^{8}}\left[4m_{c}^{3}m_{s}-2M^{2}m_{c}(m_{c}+3m_{s})\right.$
$\displaystyle\left.+2M^{4}\right]e^{-m_{c}^{2}/M^{2}}+\langle\frac{\alpha_{s}G^{2}}{\pi}\rangle\langle\overline{s}g\sigma
Gs\rangle\frac{\pi^{2}m_{c}}{108M^{12}}$
$\displaystyle\times\left[9M^{6}-60m_{c}^{5}m_{s}+10m_{c}^{3}M^{2}(3m_{c}+7m_{s})\right.$
$\displaystyle\left.-4M^{4}m_{c}(9m_{c}+4m_{s})\right]e^{-m_{c}^{2}/M^{2}}.$
(54)

Numerical computations lead to the result

$|g_{2}|=(9.8\pm 1.2)\times 10^{-1}~{}\mathrm{GeV}^{-1}.$ (55)

The partial width of the decay $M\rightarrow D_{s}^{\ast+}\rho^{+}$ is
determined by the formula

$\Gamma_{2}\left[M\rightarrow
D_{s}^{\ast+}\rho^{+}\right]=\frac{|g_{2}|^{2}m_{\rho}^{2}}{8\pi}\widehat{\lambda}\left(3+\frac{2\widehat{\lambda}^{2}}{m_{\rho}^{2}}\right),$
(56)

where $\widehat{\lambda}=\lambda(m,m_{D_{s}^{\ast}},m_{\rho})$. Then, we get

$\Gamma_{2}\left[M\rightarrow D_{s}^{\ast+}\rho^{+}\right]=(15\pm
4)~{}\mathrm{MeV}.$ (57)

The decay $M\rightarrow D^{\ast+}K^{\ast+}$ can be considered in a similar
way. Omitting details, we write down predictions for the strong coupling
$g_{3}$ and partial width of the process $M\rightarrow D^{\ast+}K^{\ast+}$:

$g_{3}=(1.5\pm 0.3)~{}\mathrm{GeV}^{-1},$ (58)

and

$\Gamma_{3}\left[M\rightarrow D^{\ast+}K^{\ast+}\right]=(37\pm
10)~{}\mathrm{MeV}.$ (59)

Information gained in this section allows us to compare couplings of the
current $J(x)$ to different two-meson states. It is seen that $g_{i}$
corresponding to vertices $MD^{\ast+}K^{\ast+}$, $MD_{s}^{\ast+}\rho^{+}$ and
$MD_{s}^{+}\pi^{+}$ obey the inequalities $g_{3}>|g_{2}|>g_{1}$. Hence, $J(x)$
describes mainly the hadronic molecule $D^{\ast+}K^{\ast+}$, as it has been
asserted in the Sec. II. Differences in the widths of relevant decays are
connected not only with $g_{i}$, but generated also by $\lambda$\- factors and
parameters of produced mesons.

Using results obtained for the partial widths of the decays considered in this
section, we can estimate the full width of the molecule $M$

$\Gamma=(123\pm 25)~{}\mathrm{MeV},$ (60)

which should be confronted with the experimental data Eq. (1). It is seen
that, although $\Gamma$ and $\Gamma_{\mathrm{exp}}$ do not coincide, they are
comparable with each other provided one takes into account errors of the
measurements and theoretical analyses.

##  IV Conclusions

In the present article, we have investigated features of the hadronic molecule
$M=D^{\ast+}K^{\ast+}$, and calculated its mass and width. The mass of $M$ has
been computed using QCD two-point sum rule method. Prediction obtained for the
mass $m=(2924~{}\pm 107)~{}\mathrm{MeV}$ is in nice agreement with the LHCb
datum for the mass of the resonance $T_{cs0}^{a++}$. It also does not differ
considerably from the mass $(2917~{}\pm 135)~{}\mathrm{MeV}$ of the molecule
$D_{s}^{\ast+}\rho^{+}$ suggested to model $T_{cs0}^{a++}$ in our paper
Agaev:2022duz  .

We have evaluated the width of the molecule $M$ by calculating partial widths
of the decay channels $M\rightarrow D_{s}^{+}\pi^{+}$, $M\rightarrow
D_{s}^{\ast+}\rho^{+}$, and $M\rightarrow D^{\ast+}K^{\ast+}$. To this end, we
have found strong couplings of the particles at vertices $MD_{s}^{+}\pi^{+}$,
$MD_{s}^{\ast+}\rho^{+}$, and $MD^{\ast+}K^{\ast+}$ by means of QCD light-cone
sum rule approach and soft-meson approximation. The strong coupling of the
hadronic molecule $M$ with final-state mesons is large in the case of mesons
$D^{\ast+}$ and $K^{\ast+}$, which is understandable because $M$ is composed
of these particles. Nevertheless, the partial width of the decay $M\rightarrow
D^{\ast+}K^{\ast+}$ is less than that of the decay $M\rightarrow
D_{s}^{+}\pi^{+}$. The reason is that the kinematical factor $\lambda$ in the
expression of the decay widths in Eqs. (40) and (56) gets its largest value in
the process $M\rightarrow D_{s}^{+}\pi^{+}$. As a result, the dominant channel
of $M$ is the decay to a pair of the mesons $D_{s}^{+}$ and $\pi^{+}$, which
was actually observed by the LHCb collaboration.

The final result $\Gamma=(123\pm 25)~{}\mathrm{MeV}$ for the full width of the
molecule $M$ is compatible with the $\Gamma_{\mathrm{exp}}$ within the
existing experimental and theoretical errors. The estimate for $\Gamma$ may be
further improved by taking into account another decay channels of the hadronic
molecule $M$. It will be interesting also to compare the $\Gamma$ and
$\Gamma_{\mathrm{exp}}$ with predictions for the full width of the resonance
$T_{cs0}^{a++}$ obtained in the context of alternative methods. Nevertheless,
based on our present results, we may consider the hadronic molecule
$M=D^{\ast+}K^{\ast+}$ as a possible candidate to the doubly charged resonance
$T_{cs0}^{a++}$.

The isoscalar partner of $T_{cs0}^{a++}$, namely the second resonance
$T_{cs0}^{a0}$ with quark content $cd\overline{s}\overline{u}$ may be modeled
as a linear superposition of hadronic molecules $D^{\ast 0}K^{\ast 0}$ and
$D_{s}^{\ast+}\rho^{-}$. The dominant decay channel of this state is the
process $T_{cs0}^{a0}\rightarrow D_{s}^{+}\pi^{-}$. The modes
$T_{cs0}^{a0}\rightarrow D^{0}K^{0}$, $D_{s}^{\ast+}\rho^{-}$ and $D^{\ast
0}K^{\ast 0}$ are other possible decay channels of $T_{cs0}^{a0}$.
Experimentally measured mass and width differences between $T_{cs0}^{a++}$ and
$T_{cs0}^{a0}$ are equal to $\Delta m\approx 28~{}\mathrm{MeV}$ and
$\Delta\Gamma\approx 15~{}\mathrm{MeV}$, respectively. To be accepted as a
reliable model for $T_{cs0}^{a0}$ the molecule picture should be successfully
confronted with the available data.

Another problem to be addressed here, is similarities and differences of the
resonances $T_{cs0}^{a0/++}$ and $X_{0(1)}$. It has been noted in Sec. I that
$X_{0(1)}$ are exotic mesons composed of quarks $ud\overline{s}\overline{c}$.
The scalar structures $T_{cs0}^{a0}$ and $X_{0}$ differ from each other by
quark-exchanges $\overline{u}\leftrightarrow u$ and
$c\leftrightarrow\overline{c}$. They are neutral particles with masses
$2892~{}\mathrm{MeV}$ and $2866~{}\mathrm{MeV}$, respectively. It is seen,
that a mass gap between these two structures is small. In Ref. Agaev:2020nrc ,
we modeled $X_{0}$ as a hadronic molecule $\overline{D}^{\ast 0}K^{\ast 0}$
and found its mass equal to $2868~{}\mathrm{MeV}$, which is in a very nice
agreement with the LHCb datum. The $\overline{D}^{\ast 0}K^{\ast 0}$ and a
component $D^{\ast 0}K^{\ast 0}$ of the molecule model for $T_{cs0}^{a0}$,
have almost identical structures, therefore one expects the mass of the
isoscalar state $T_{cs0}^{a0}$ will be be consistent with experiments.

The mass difference $\sim 55~{}\mathrm{MeV}$ between $T_{cs0}^{a++}$ and
$X_{0}$ is larger than the gap in the previous case which is connected
presumably with doubly charged nature of $T_{cs0}^{a++}$. We explored the
vector resonance $X_{1}$ in Ref. Agaev:2021knl  as a diquark-antidiquark state
$[ud][\overline{c}\overline{s}]$ and achieved reasonable agreements with the
LHCb data. In general, hadronic molecules composed of two mesons may be used
to model vector particles as well. Because $T_{cs0}^{a++}$ has the spin-parity
$J^{\mathrm{P}}=0^{+}$, in this article, we have studied only scalar particle
$D^{\ast+}K^{\ast+}$. Molecules with the same quark content but different
spin-parities are yet hypothetical structures. They are interesting objects
for theoretical researches as well, because may be discovered soon in various
exclusive processes. Properties of the isoscalar resonance $T_{cs0}^{a0}$, as
well as counterparts of $T_{cs0}^{a0/++}$ with different spin-parities are
issues for future investigations.

## References

  * (1) [LHCb], arXiv:2212.02716 [hep-ex].
  * (2) [LHCb], arXiv:2212.02717 [hep-ex].
  * (3) R. Aaij et al. [LHCb], Phys. Rev. Lett. 125, 242001 (2020).
  * (4) R. Aaij et al. [LHCb], Phys. Rev. D 102, 112003 (2020).
  * (5) V. M. Abazov et al. [D0], Phys. Rev. Lett. 117, 022003 (2017).
  * (6) V. M. Abazov et al. [D0], Phys. Rev. D 97, 092004 (2018).
  * (7) S. S. Agaev, K. Azizi and H. Sundu, Phys. Rev. D 93, 094006 (2016).
  * (8) W. Chen, H. X. Chen, X. Liu, T. G. Steele and S. L. Zhu, Phys. Rev. Lett. 117, 022002 (2016).
  * (9) S. S. Agaev, K. Azizi and H. Sundu, Eur. Phys. J. C 78, 141 (2018).
  * (10) S. S. Agaev, K. Azizi and H. Sundu, Phys. Lett. B 820, 136530 (2021).
  * (11) M. Karliner and J. L. Rosner, Phys. Rev. D 102, 094016 (2020).
  * (12) Z. G. Wang, Int. J. Mod. Phys. A 35, 2050187 (2020).
  * (13) H. X. Chen, W. Chen, R. R. Dong and N. Su, Chin. Phys. Lett. 37, 101201 (2020).
  * (14) M. Z. Liu, J. J. Xie and L. S. Geng, Phys. Rev. D 102, 091502 (2020).
  * (15) X. H. Liu, M. J. Yan, H. W. Ke, G. Li and J. J. Xie, arXiv:2008.07190 [hep-ph].
  * (16) R. Molina and E. Oset, Phys. Lett. B 811, 135870 (2020).
  * (17) M. W. Hu, X. Y. Lao, P. Ling and Q. Wang, Chin. Phys. C 45, 021003 (2021).
  * (18) X. G. He, W. Wang and R. Zhu, Eur. Phys. J. C 80, 1026 (2020).
  * (19) Q. F. Lu, D. Y. Chen and Y. B. Dong, Phys. Rev. D 102, 074021 (2020).
  * (20) J. R. Zhang, Phys. Rev. D 103, 054019 (2021).
  * (21) Y. Huang, J. X. Lu, J. J. Xie and L. S. Geng, Eur. Phys. J. C 80, 973 (2020).
  * (22) Y. Xue, X. Jin, H. Huang and J. Ping, Phys. Rev. D 103, 054010 (2021).
  * (23) G. Yang, J. Ping and J. Segovia, Phys. Rev. D 103, 074011 (2021).
  * (24) T. W. Wu, Y. W. Pan, M. Z. Liu, S. Q. Luo, L. S. Geng, and X. Liu, Phys. Rev. D 105, L031505 (2022).
  * (25) L. M. Abreu, Phys. Rev. D 103, 036013 (2021).
  * (26) G. J. Wang, L. Meng, L. Y. Xiao, M. Oka and S. L. Zhu, Eur. Phys. J. C 81, 188 (2021).
  * (27) C. J. Xiao, D. Y. Chen, Y. B. Dong and G. W. Meng, Phys. Rev. D 103, 034004 (2021).
  * (28) X. K. Dong and B. S. Zou, Eur. Phys. J. A 57, 139 (2021).
  * (29) T. J. Burns and E. S. Swanson, Phys. Rev. D 103, 014004 (2021).
  * (30) A. E. Bondar and A. I. Milstein, JHEP 12, 015 (2020).
  * (31) Y. K. Chen, J. J. Han, Q. F. Lü, J. P. Wang and F. S. Yu, Eur. Phys. J. C 81, 71 (2021).
  * (32) R. M. Albuquerque, S. Narison, D. Rabetiarivony and G. Randriamanatrika, Nucl. Phys. A 1007, 122113 (2021).
  * (33) B. Wang, and S. L. Zhu, Eur. Phys. J. C 82, 419 (2022).
  * (34) S. S. Agaev, K. Azizi and H. Sundu, J. Phys. G 48, 085012 (2021).
  * (35) S. S. Agaev, K. Azizi and H. Sundu, Nucl. Phys. A 1011, 122202 (2021).
  * (36) S. S. Agaev, K. Azizi and H. Sundu, Phys. Rev. D 106, 014019 (2022).
  * (37) H. Sundu, S. S. Agaev, and K. Azizi, Eur. Phys. J. C 83, 198 (2023).
  * (38) S. S. Agaev, K. Azizi and H. Sundu, J. Phys. G to appear, arXiv:2207.02648 [hep-ph].
  * (39) R. Chen, and Q.Huang, arXiv:2208.10196 [hep-ph].
  * (40) Y. H. Ge, X. H. Liu,and H. W. Ke, Eur. Phys. J. C 82, 955 (2022).
  * (41) J. Wei, Y. H. Wang, C. S. An, and C. R. Deng, Phys. Rev. D 106, 096023 (2022).
  * (42) R. Molina, and E. Oset, Phys. Rev. D 107, 056015 (2023).
  * (43) F. X. Liu, R. H. Ni, X. H. Zhong, and Q. Zhao, arXiv:2211.01711 [hep-ph].
  * (44) Z. L. Yue, C. J. Xiao, and D. Y. Chen, arXiv:2212.03018 [hep-ph].
  * (45) Q. Qin, J. L. Qui, and F. S. Yu, Eur. Phys. J. C 83, 227 (2023).
  * (46) M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B 147, 385 (1979).
  * (47) M. A. Shifman, A. I. Vainshtein and V. I. Zakharov, Nucl. Phys. B 147, 448 (1979).
  * (48) I. I. Balitsky, V. M. Braun and A. V. Kolesnichenko, Nucl. Phys. B 312, 509 (1989).
  * (49) B. L. Ioffe and A. V. Smilga, Nucl. Phys. B 232, 109 (1984).
  * (50) V. M. Belyaev, V. M. Braun, A. Khodjamirian and R. Ruckl, Phys. Rev. D 51, 6177 (1995).
  * (51) Q. N. Wang, W. Chen, and H. X. Chen, Chin. Phys. C 45, 093102 (2021).
  * (52) H. X. Chen, Y. X. Yan, W. Chen, Phys. Rev. D 106, 094019 (2022).
  * (53) Q. Xin, and Z. G. Wang, Eur. Phys. J. A 58, 118 (2022).
  * (54) S. S. Agaev, K. Azizi and H. Sundu, Turk. J. Phys. 44, 95 (2020).
  * (55) R. L. Workman et al. [Particle Data Group], Prog. Theor. Exp. Phys. 2022, 083C01 (2022).
  * (56) S. S. Agaev, K. Azizi and H. Sundu, Phys. Rev. D 93, 074002 (2016).
  * (57) S. S. Agaev, K. Azizi and H. Sundu, Phys. Rev. D 106, 014025 (2022).
  * (58) S. S. Agaev, K. Azizi and H. Sundu, Phys. Rev. D 93, 114036 (2016).
  * (59) V. Lubicz, A. Melis and S. Simula, PoS LATTICE 2016, 291 (2017).