\section{Introduction}

In the past few decades, especially for the discovery of the resonance $a_0(980)$~\cite{Ammar:1968zur}, the properties of scalar mesons still remain the controversial, which become hot topics in hadron physics. Theoretically, they have been considered as tetraquark states~\cite{Jaffe:1976ig, Alford:2000mm, Humanic:2022hpq, Brito:2004tv, Klempt:2007cp} or two-meson bound states~\cite{Weinstein:1982gc, Branz:2007xp, Dai:2012kf, Dai:2014lza, Sekihara:2014qxa}. In the literature, $a_0(980)$ is usually taken as a scalar meson with $q\bar q$ configuration based on the constituent quark model. In the present paper, we shall adopt this $q\bar q$ configuration for $a_0(980)$ and make a detailed study on its twist-2 light-cone distribution amplitudes, which can be tested when more accurate data have been achieved in the near future.

The semileptonic decay of the charmed meson provides simpler decay mechanism and final-state interactions, and it is an ideal platform for studying the meson's properties. The BESIII collaboration observed the charmless hadronic decay processes involving $a_0(980)$-meson, and reported the results of $D^0\to a_0 (980)^- e^+\nu_e$ and $D^+\to a_0 (980)^0e^+ \nu_e$ with the significance up to $6.4\sigma$ and $2.9\sigma$, respectively~\cite{BESIII:2018sjg}. The branching ratio of $D\to a_0(980)\ell\bar \nu_\ell$ combing with $a_0(980)$-meson decays $a_0(980)\to \eta \pi$ are ${\mathcal B}(D^0\to a_0(980)^- e^+\nu_e)\times{\mathcal B}(a_0(980)^-\to \eta \pi^-)=(1. 33^{+0.33}_{-0.29}\pm0.09)\times10^{-4}$ and ${\cal B}(D^+\to a_0(980)^0 e^+ \nu_e)\times {\cal B} (a_0(980)^0\to \eta \pi^0) = (1.66^{+0.81}_{-0.66} \pm0.11) \times 10^{-4}$, where the first and second errors are statistical error and systematic error, respectively. Recently, the BESIII collaboration presented the branching ratio of the decay $D_s^+ \to f_0(980)(\to \pi^0 \pi^0) e^+\nu_e$ for the first time~\cite{BESIII:2022dxg}, by using the $6.32~{\rm fb}^{-1}$ $e^+e^-$ collision data recorded with the center-of-mass energies between 4.178 GeV to 4.226 GeV. These experimental values will help us to study the production mechanism of scalar mesons through the hadron decay process. The $D\to a_0(980)\ell\nu_\ell$ semileptonic decays also provide ideal platforms for deep learning of the $a_0(980)$-meson and its related preferences. And the important parameters for the semileptonic decays $D\to a_0(980)\ell\bar \nu_\ell$, many theoretical methods have been employed to calculate the transition form factors (TFFs), such as the covariant confined quark model (CCQM)~\cite{Soni:2020sgn}, the light-cone sum rules (LCSR)~\cite{Cheng:2017fkw, Huang:2021owr} and the anti-de Sitter with QCD potential AdS/QCD~\cite{Momeni:2022gqb}.

The LCSR is an effective method in dealing with the heavy to light decays. It applies a partial resummation of the operator product expansion to all orders and reorganizes the operator production expansion in terms of the twists of relevant operators, which then results in a series over the light-meson's light-cone distribution amplitudes (LCDAs) of increasing twists. The LCSR is applicable in a wide $q^2$-region and could be extended to all $q^2$-regions via proper extrapolations. In this paper, we will adopt the LCSR to calculate the $D\to a_0(980)$ TFFs. One of the key quantities that characterize the TFFs is the twist-2 LCDA, which describes the dominant momentum fraction distribution for each part of a meson. In general, the $a_0(980)$-meson twist-2 LCDA $\phi_{2;a_0}(x,\mu)$ can be expanded as a Gegenbauer polynomial series~\cite{Cheng:2005nb}
\begin{eqnarray}
\phi_{2;a_0}(x,\mu) = 6x\bar x \bigg[a_{2;a_0}^0(\mu)  + \sum\limits_{n=1}^\infty a_{2;a_0}^n(\mu) C_n^{3/2}(\xi) \bigg],
\end{eqnarray}
where $\bar x = (1-x)$ and $\xi = (2x-1)$. The even Gegenbauer coefficients are highly suppressed, which exactly equal to zero under the approximation that $m_1\simeq m_2$ ($m_{1,2}$ are masses of two constituent quarks), and the LCDA of the scalar meson is then dominated by the odd Gegenabuer moments. In contrast, the odd Gegenbauer moments vanish for $\pi$ and $\rho$ mesons. Thus, it tends to the antisymmetric form under $u\to (1-u)$ transition in the SU$_f$(3) limit. Currently, the $a_0(980)$-meson twist-2 LCDA is mainly coming from QCDSR by Cheng, Chua and Yang~\cite{Cheng:2005nb}, which gives the first two order Gegenbauer moments. To have a deeper insight into the $a_0(980)$-meson twist-2 LCDA and to improve the accuracy of the processes involving $a_0(980)$-meson, it is helpful to known its higher order Gegenbauer moments, which can be achieved by establishing a light-cone harmonic oscillator (LCHO) model with the help of the Brodsky-Huang-Lepage (BHL) prescription. The LCHO model parameters can be fixed by fitting them to the known moments $\langle\xi^n_{2;a_0}\rangle|_\mu$ at a specified scale $\mu$ with the least squares method.

For the purpose, we need to calculate the $a_0(980)$-meson twist-2 LCDA moments $\langle\xi^n_{2;a_0}\rangle|_\mu$. As one of the effective way to calculate each order moments for light or heavy meson twist-2 or twist-3 LCDA is the QCDSR under the framework of the background field theory (BFTSR). The basic assumption of QCDSR is the introduction of vacuum condensates. In addition, QCD background field method systematically describes these vacuum condensates from the perspective of field theory~\cite{Novikov:1983gd, Hubschmid:1982pa, Govaerts:1983ka, Govaerts:1984bk, Ambjorn:1982bp, Ambjorn:1982en, Reinders:1984sr, Elias:1987ac, Huang:1986wm, Huang:1989gv}. Quark and gluon fields are composed of background fields and quantum fluctuations around them under the framework of BFTSR. The vacuum condensates and quantum fluctuations of the background field represent non perturbative effects and computable perturbative effects respectively. These vacuum condensates provide a clear physical image for the internal structure of the bound state, which makes the calculation easier. At present, the BFTSR has been applied to calculate the LCDAs of pseudoscalar and vector/axial vector mesons~\cite{Fu:2016yzx,Hu:2021zmy, Hu:2021lkl, Zhong:2014jla, Fu:2018vap, Zhong:2016kuv, Huang:2004tp, Zhong:2011rg, Zhong:2022lmn}. In this paper, we will study the scalar $a_0(980)$-meson twist-2 LCDA by using BFTSR firstly.

The remaining parts of this paper are organized as follows. In Sec.~\ref{Sec:II}, we calculate the $a_0(980)$-meson twist-2 LCDA moments and introduce the LCHO model. A new improved model is proposed and the model parameters shall be obtained by fitting moment with the least square method. Section~\ref{Sec:III} gives the numerical results and discussions. Section~\ref{Sec:IV} is for a brief summary.

\section{Theoretical framework}\label{Sec:II}

The $a_0(980)$-meson have three types of states, one is $a_0(980)^0$-meson with $(\bar uu-\bar dd)/\sqrt 2$ component, the other is $a_0(980)^-$-meson with $(\bar ud)$ component, and the third one is $a_0(980)^+$-meson with $(\bar du)$. Based on the basic procedure of QCD sum rules, one can adopt the following two-point correlator to derive the sum rules for $a_0(980)$-meson twist-2 LCDA moments $\langle\xi^n_{2;a_0}\rangle|_\mu$, which can be read off,
\begin{eqnarray}
\Pi^{(n,0)}_{2;a_0} (z,q) &=& i \int d^4x e^{iq\cdot x} \langle 0| T \{ J_n^V(x), J^{S,\dag}_0(0) \} |0\rangle
\nonumber\\
&=& (z\cdot p)^{n+1} I^{(n,0)}_{2;a_0} (q^2),
\label{Eq:correlator}
\end{eqnarray}
where $z^2 = 0$ and $n$ take odd numbers. The currents are taken as $J^V_n(x) = \bar q_1(x) \DS z (i z\cdot \tensor{D})^n q_2(x)$ and $J^S_0(0) = \bar q_1(0)q_2(0)$ with the covariant derivative $(iz \cdot \tensor D)^n=(iz\cdot \overrightarrow D - iz\cdot \overleftarrow D)^n$, which are mainly come from definitions of the scalar $a_0(980)$-meson twist-2 LCDA $\phi_{2;a_0}(x,\mu)$ and twist-3 LCDA $\phi_{3,a_0}^p (x,\mu)$, which can be written as~\cite{Cheng:2005nb}
\begin{align}
&\hspace{-0.2cm}\langle 0|\bar q_1(z)\gamma _\mu q_2(-z)|a_0(p)\rangle
\nonumber\\
&\qquad\qquad = p_\mu f_{a_0} \int^1_0 d xe^{i(2u-1)(p\cdot z)}\phi_{2;a_0}(x,\mu),
\label{Eq:a0}
\\
&\hspace{-0.2cm}\langle 0|\bar q_1(z) q_2(-z)|a_0(p)\rangle
\nonumber\\
&\qquad\qquad = m_{a_0}f_{a_0} \int^1_0 dx e^{i(2u-1)(p\cdot z)}\phi_{3,a_0}^p(x,\mu),
\label{Eq:phi3s}
\end{align}
 From one side, one can apply the OPE for the correlator, e.g. Eq.~\eqref{Eq:correlator} in deep Euclidean region $q^2 \ll 0$. Then, the correlator can be expanded into three terms including the quark propagators and the vertex operators,
\begin{align}
\Pi^{(n,0)}_{2;a_0}(z,q)&= i\int d^4x e^{iq\cdot x}\Big\{
\nonumber\\
&-{\rm Tr}\langle 0|S^{q_1}_{F}(0,x)\DS z (iz\cdot\tensor D)^n S^{q_2}_F(x,0)|0\rangle
\nonumber\\
& +{\rm Tr}\langle 0|\bar q_1(x)q_1(0)\DS z(iz\cdot\tensor{D})^n S^{q_2}_F(x,0)|0\rangle
\nonumber\\
& +{\rm Tr}\langle 0|S^{q_1}_F(0,x)\DS z(iz\cdot\tensor{D})^n \bar q_2(x)q_2(0)|0\rangle
\nonumber\\
& +\cdots \Big\}
\label{pp}
\end{align}
In dealing with Lorentz invariant scalar function $\Pi_{2;a_0}^{(n,0)}(z, q^2)$, the vacuum matrix element should be used, which can be found in our previous work~\cite{Zhong:2014jla}. Meanwhlie, we take $q_1=q_2=q$ stand for the light $u$,$d$-quark. On the other hand, one can insert a complete set of $a_0(980)$-meson intermediated hadronic states with the same $J^P$ quantum number into the correlator and obtain the hadronic expression
\begin{align}
\textrm{Im} I^{(n,0)}_{2;a_0,{\rm had}}(q^2)&=\pi \delta (q^2 - m^{2}_{a_0}) f^{2}_{a_0} \langle\xi^n_{2;a_0}\rangle|_\mu\langle \xi^{p;0}_{3;a_0}\rangle|_\mu
\nonumber\\
&+\frac{3m_q}{4\pi^2(n+2)}\theta(q^2- s_{a_0}).
\label{hadim}
\end{align}
In which the definitions of moments for $a_0(980)$-meson twist-2, 3 LCDAs are used, which have the following formula
\begin{align}
& \hspace{-0.5cm}\langle 0 |\bar q_1 (0) \DS z (iz \cdot\tensor D)^n q_2(0)|a_0(p)\rangle = (z\cdot p)^{n+1}f_{a_0}\langle\xi^n_{2;a_0}\rangle|_\mu,
\label{DA2}\\
& \hspace{-0.5cm} \langle 0|\bar q_1(0) q_2(0)|a_0(p)\rangle =m_{a_0}f_{a_0}\langle\xi^{p,0}_{3;a_0}\rangle|_\mu.
 \label{DA3}
 \end{align}
where $m_q=m_{q_1}=m_{q_2}$ in Eqs.~\eqref{hadim}-\eqref{DA3} stand for the mass of light quark and the $m_{a_0}$ is the mass of $a_0(980)$-meson, $f_{a_0}$ is the decay constant. The slight difference between mass of $u$-quark and $d$-quark is ignored in this paper. Meanwhile, the $s_{a_0}$ stands for the continuum threshold. After matching the hadronic and OPE sides through the dispersion relation and performing the Borel transformation, one can get the sum rule for $\langle\xi_{2;a_0}^n\rangle|_\mu \langle\xi_{3;a_0}^{p;0}\rangle|_\mu$,
\begin{align}
\langle\xi^n_{2;a_0}\rangle|_\mu\langle \xi^{p;0}_{3;a_0}\rangle|_\mu
&= I_1^{\rm pert.}(n,s_0,M^2)
\nonumber\\
&+ \langle\bar qq \rangle I_1^{\langle\bar qq \rangle}(n,M^2)
\nonumber\\
&+ \langle g_s\bar q q \rangle^2 I_1^{\langle g_s\bar q q \rangle^2}(n,M^2)
\nonumber\\
& + \langle g_s^2\bar qq\rangle^2 I_1^{\langle g_s^2\bar qq\rangle^2}(n,M^2,\mu)
\nonumber\\
&+ \langle\alpha_s G^2\rangle I_1^{\langle\alpha_s G^2\rangle}(n,M^2,\mu)
\nonumber\\
& + \langle g_s^3fG^3\rangle I_1^{\langle g_s^3fG^3\rangle}(n,M^2,\mu)
\nonumber\\
&+ \langle{g_s}\bar q \sigma TGq\rangle I_1^{\langle g_s\bar q \sigma TGq\rangle}(n,M^2).
\label{Eq:SRxi2nxi30}
\end{align}
In which, we separate the sum rules into seven different terms. The detailed expressions for $I_1^{i}(n, s_0, M^2, \mu)$ with $i=({\rm pert.}, {\langle\bar qq \rangle}, {\langle g_s\bar q q \rangle^2}, {\langle\alpha_s G^2\rangle}, {\langle g_s^3fG^3\rangle}, {\langle g_s\bar q \sigma TGq\rangle})$ accompany by the six different vacuum condensates and perturbative terms are read off:
\begin{widetext}
\begin{eqnarray}
&&\hspace{-0.5cm}
  I_1^{\rm pert.} (n,M^2) = \frac{m_q}{m_{a_0}f_{a_0}^2} e^{m_{a_0}^2/M^2}\frac{3}{8\pi^2(n+2)} \left(1-e^{-s_0/M^2}\right),
\label{Eq:B1pert}
\\
&&\hspace{-0.5cm}
  I_1^{\langle\bar qq \rangle}(n,M^2) = \frac{2}{m_{a_0}f_{a_0}^2} e^{m_{a_0}^2/M^2},
\label{Eq:B2qq}
\\
&&\hspace{-0.5cm}
  I_1^{\langle g_s\bar q q \rangle^2}(n,M^2) = \frac{m_q}{m_{a_0}f_{a_0}^2} e^{m_{a_0}^2/M^2}\frac{4(n+3)}{81M^4},
\label{Eq:B3gsqq2}
\\
&&\hspace{-0.5cm}
   I_1^{\langle g_s^2\bar qq\rangle^2}(n,M^2,\mu) = -\frac{2m_q}{m_{a_0}f_{a_0}^2} e^{m_{a_0}^2/M^2}\frac{2+\kappa^2}{3888\pi^2 M^4} \bigg\{\delta^{0n}\Big[-24\bigg(-\ln\frac{M^2}{\mu^2}\bigg)-\,148\Big]+\delta^{1n}\bigg[128\bigg(-\,\ln\frac{M^2}{\mu^2}\bigg)-\,692\bigg]
  \hspace{0.7cm}
\nonumber\\
  &&\qquad+\theta(n-1)\bigg[8(6n^2+34n)(-\!\ln\frac{M^2}{\mu^2})\!+\!4n\tilde\psi(n)-2(6n^2+96n+212)\bigg]+\theta(n-2)\bigg[8(33n^2-17n)
  \bigg(\!-\!\ln\frac{M^2}{\mu^2}\bigg)
\nonumber\\
  &&\qquad\!-\!2(6n^2\!+71n)\tilde\psi(n)\!-\!\frac{1}{n(n-1)}(231n^4+520n^3-1101n^2+230n)\bigg]+\theta(n-3)
  \bigg[(74n-144n^2)\tilde\psi(n)-\frac{1}{n-1}
\nonumber\\
    &&\qquad\times(169n^3-348n^2+245n+60)\bigg]+4(n+5)\bigg\},
\label{Eq:B4gs2qq2}
\\
&&\hspace{-0.5cm}
    I_1^{\langle\alpha_s G^2\rangle}(n,M^2,\mu) = -\!\frac{2m_q}{m_{a_0}f_{a_0}^2}e^{m_{a_0}^2/M^2}\frac1{48\pi M^2 }\bigg\{12n\bigg(-\!\!\ln\frac{M^2}{\mu^2}\bigg)\!\!-\!\!6(n+2)\!\!+\!\!\theta(n-1)\bigg[4n\bigg(-\!\ln\frac{M^2}{\mu^2}\bigg)\!+\!3\tilde\psi(n)-\frac6n\bigg]
\nonumber\\
    &&\qquad+\theta(n-2)\bigg[-(8n+3)\,\tilde\psi(n)\,-\, 2(2n+1)+ \frac6n\bigg]\bigg\},
\label{Eq:B5alphasG2}
\\
&&\hspace{-0.5cm}
    I_1^{\langle g_s^3fG^3\rangle}(n,M^2,\mu) = -~\frac{2m_q}{m_{a_0}f_{a_0}^2} e^{m_{a_0}^2/M^2}\frac{1}{384\pi^2M^4} \bigg\{\delta^{1n}\bigg[-\,24\bigg(\,-\,\,\ln\frac{M^2}{\mu^2}\bigg)\,+\,84\bigg]~+~\,\theta(n-1)\bigg[-~4n(3n-5)
\nonumber\\
    &&\qquad\times\bigg(-\ln\frac{M^2}{\mu^2}\bigg) +2(2n^2+5n-13)\bigg]\!+\!\theta(n-2)\bigg[-24n^2\bigg(-\!\ln\frac{M^2}{\mu^2}\bigg)+2n(n-4)\tilde\psi(n)+17n^2+55n+12\bigg]
\nonumber\\
    &&\qquad +\theta(n-3)\bigg[2n(n-4)\tilde\psi(n)+\frac{1}{n-1}(19n^3-32n^2+7n+6)\bigg] \!\bigg\},
\label{Eq:B5gs3fG3}
\\
&&\hspace{-0.5cm}
    I_1^{\langle g_s\bar q \sigma TGq\rangle}(n,M^2) =- \frac{1}{m_{a_0}f_{a_0}^2} e^{m_{a_0}^2/M^2} \frac{4n}{3M^2},
\end{eqnarray}
\end{widetext}
with $\tilde\psi(n)=\psi(\frac{n+1}2) - \psi(\frac n2) +(-1)^n\ln4$. It is worth to point out that the detailed terms ${\cal B}_{M^2} I_{ijk(n)}$ are listed in the Appendix~\ref{Appendix:3}. For a more thorough considering the sum rule for $\langle\xi^n_{2;a_0}\rangle|_\mu$, it might be convenient to calculate $\langle \xi^{p;0}_{3;a_0}\rangle|_\mu$. To achieve this target, one can use the correlator.
\begin{align}
\Pi^{(0,0)}_{3;a_0} (z,q) = i \int d^4x e^{iq\cdot x} \langle 0| T \{ J^S_0(x), J^{S\dag}_0(0)\} |0\rangle,
\end{align}
for the $a_0(980)$-meson twist-3 LCDA $0_{\rm th}$ moment. Followed by the basic procedure of the QCDSR, the expression of the $0_{\rm th}$ moment of scalar meson $a_0(980)$-meson twist-3 LCDA can be obtained,
\begin{align}
\langle\xi^{p;0}_{3;a_0}\rangle|_\mu
&= I_2^{\rm pert.}(n,s_0,M^2)
\nonumber\\
&+ \langle\bar qq \rangle I_2^{\langle\bar qq \rangle}(M^2)
\nonumber\\
&+ \langle g_s\bar q q \rangle^2 I_2^{\langle g_s\bar q q \rangle^2}(M^2)
\nonumber\\
& + \langle g_s^2\bar qq\rangle^2 I_2^{\langle g_s^2\bar qq\rangle^2}(M^2,\mu)
\nonumber\\
&+ \langle\alpha_s G^2\rangle I_2^{\langle\alpha_s G^2\rangle}(M^2)
\nonumber\\
&+ \langle g_s \bar q \sigma TGq\rangle I_2^{\langle g_s\bar q \sigma TGq\rangle}(M^2).
\label{Eq:xi30xi30}
\end{align}
In which, the light quark is taken as $q=(u,d)$ and corresponding vacuum condensates are given in the next Section. The detailed expressions for the subitems in Eq.~\eqref{Eq:xi30xi30} are present as follows,
\begin{eqnarray}
&&\hspace{-0.5cm}
    I_2^{\rm pert.} (s_0,M^2) = \frac{3e^{m_{a_0}^2/M^2}M^2}{8\pi^2 m_{a_0}^2 f_{a_0}^2}  \bigg[M^2\!-\!(M^2\!+\!s_0)e^{- s_0/M^2}\bigg],
\nonumber\\
\label{Eq:B1pert}
\\
&&\hspace{-0.5cm}
    I_2^{\langle\bar qq \rangle}(M^2) = \frac{3m_qe^{m_{a_0}^2/M^2}}{m_{a_0}^2 f_{a_0}^2},
\label{Eq:B2qq}
\\
&&\hspace{-0.5cm}
   I_2^{\langle g_s\bar q q \rangle^2}(M^2) = - \frac{8 e^{m_{a_0}^2/M^2}}{27 m_{a_0}^2 f_{a_0}^2 M^2},
\label{Eq:B3gsqq2}
\\
&&\hspace{-0.5cm}
   I_2^{\langle g_s^2\bar qq\rangle^2}\!(M^2,\mu) \!=\! \frac{e^{m_{a_0}^2 / M^2}\!(2\! +\!\kappa^2)} {486\pi^2 m_{a_0} f_{a_0}^2 M^2}\!\bigg[35\!-\!6\left(\!-\!\ln \frac{M^2}{\mu^2}\!\right)\!\bigg],
\nonumber\\
\label{Eq:B4gs2qq2}
\\
&&\hspace{-0.5cm}
  I_2^{\langle\alpha_s G^2\rangle}(M^2) = \frac{e^{m_{a_0}^2/M^2}}{m_{a_0}^2 f_{a_0}^2} \frac1{8\pi},
\label{Eq:B5alphasG2}
\\
&&\hspace{-0.5cm}
  I_2^{\langle g_s\bar q \sigma TGq\rangle}(M^2) = \frac{e^{m_{a_0}^2/M^2}}{m_{a_0}^2 f_{a_0}^2 M^2}.
\end{eqnarray}
Since higher-order and higher-dimensional corrections are difficult to calculate completely, the moments $\langle\xi^n_{2;a_0}\rangle|_\mu$ of the sum rule (\ref{Eq:SRxi2nxi30}) cannot be normalized in the whole Borel parameter region. To get a more accurate moments $\langle\xi^n_{2;a_0}\rangle|_\mu$ for sum rule, we can use the following expression
\begin{align}
\langle\xi^n_{2;a_0}\rangle|_\mu=\frac{\langle\xi^n_{2;a_0}\rangle|_\mu\langle\xi^{p;0}_{3;a_0}\rangle|_\mu}{\sqrt{\langle\xi^{p;0}_{3;a_0}
\rangle^2|_\mu}}
\end{align}
This method can eliminate the systematic errors caused by many factors. The discussion for the pion and kaon cases can be found in our previous work~\cite{Zhong:2021epq,Zhong:2022ecl}.

Normally, conformal expansion of LCDAs Gegenbauer polynomials makes the higher-order Gegenbauer moments unreliable. To improve this situation, the LCHO is adopted to determine $a_0(980)$-meson twist-2 LCDA. Referring to the LCHO model of the pion leading-twist WF raised in Refs.~\cite{Wu:2010zc, Wu:2011gf}, it is expressed as:
\begin{align}
&\Psi_{2;a_0}(x,\textbf{k}_\bot)= \sum_{\lambda_1\lambda_2} \chi_{2;a_0}^{\lambda_1\lambda_2}(x,\textbf{k}_\bot) \Psi^R_{2;a_0}(x,\textbf{k}_\bot)
\label{WF_full}
\end{align}
where $\textbf{k}_\bot$ is transverse momentum. Furthermore, $\lambda_1$ and $\lambda_2$ are the helicities of the two constituent quark. $\chi_{2;a_0} ^{\lambda_1 \lambda_2} (x, \textbf{k}_\bot)$ stands for the spin-space WF that comes from the Wigner-Melosh rotation. The different forms of $\lambda_1\lambda_2$ are shown in Table \ref{table:chi}, which can also been seen in Refs.~\cite{Huang:1994dy, Cao:1997hw, Huang:2004su, Wu:2005kq}. The spin-space WF, e.g. $\sum_{\lambda_1 \lambda_2} \chi_{2;a_0}^{\lambda_1 \lambda_2} (x,\textbf{k}_\bot) ={\hat m_q}^2/(\textbf{k}^2_\bot + {\hat m_q}^2)^{1/2}$. Then, by combing the spatial WF $\Psi^R_{2;a_0}(x,\textbf{k}_\bot)=A_{2;a_0} \varphi_{2;a_0}(x) \exp[ -(\textbf{k}^2_\bot + {\hat m_q}^2)/(8\beta_{2;a_0}^2 x\bar x)]$, the $a_0(980)$-meson WF will be obtained. In which, the $A_{2;a_0}$, $\hat m_q$ stand for the normalization constant and light quark mass. The final $a_0(980)$-meson twist-2 LCDA can be obtained by using the relationship between the twist-2 LCDA and WF of $a_0(980)$-meson, e.g integrated over the squared transverse momentum, which can be expressed as
\begin{align}
&\phi_{2;a_0}(x,\mu)=\frac{ A_{2;a_0}\hat m_q \beta_{2;a_0}}{4\sqrt{2}\pi^{3/2}} \sqrt{x\bar x} \varphi_{2;a_0}(x)
\nonumber\\
&\quad\times \left\{ \textrm{Erf}\left[ \sqrt{\frac{{\hat m_q}^2+\mu^2 }{8\beta_{2;a_0}^2 x\bar x}} \right] - \textrm{Erf}\left[ \sqrt{\frac{{\hat m_q^2}}{8\beta_{2;a_0}^2 x\bar x}} \right] \right\},
\label{DA_model}
\end{align}
where ${\rm Erf}(x) = 2\int^x_0 e^{-t^2} dx/{\sqrt{\pi}}$ is the error function, and $\varphi_{2;a_0}(x) = (x\bar x)^{\alpha_{2;a_0}} C_1^{3/2}(2x-1)$. Based on the experience of other mesons~\cite{Huang:2013gra, Huang:2013yya, Wu:2012kw, Zhong:2015nxa, Zhang:2021wnv, Hu:2021lkl, Zhong:2018exo, Zhong:2022ecl, Zhong:2021epq}, we take the wavefunction parameter $\beta_{2;a_0}=0.5$. Whether the value of $\beta_{2;a_0}$ is accurate can be judged by goodness of fit $P_{\chi^2}$. The free parameters $\alpha_{2;a_0}$ and $A_{2;a_0}$ can be obtained by fitting the moments $\langle \xi^n_{2;a_0} \rangle|_\mu$ directly.

\begin{table}
\centering
\renewcommand\arraystretch{1.3}
\small
\caption{The expressions of the spin-space wave function $\chi_{2;a_0}^{\lambda_1 \lambda_2}(x,\textbf{k}_\perp)$ with different $\lambda_1 \lambda_2$.}\label{table:chi}
\begin{tabular}{c c c c}
\\ \hline
$\lambda_1 \lambda_2$ & $\chi_{2;a_0}^{\lambda_1 \lambda_2}(x,\textbf{k}_\perp)$ &
~~~~~~~~~~~~~$\lambda_1 \lambda_2$ & $\chi_{2;a_0}^{\lambda_1 \lambda_2}(x,\textbf{k}_\perp)$\\
\hline
$\downarrow\downarrow$ & $-\dfrac{k_x + i k_y}{\sqrt{2({\hat m_q}^2 + \textbf{k}_\perp^2)}}$ &
~~~~~~~~~~~~~$\uparrow\uparrow$   & $-\dfrac{k_x - i k_y}{\sqrt{2({\hat m_q}^2 + \textbf{k}_\perp^2)}}$ \\
$\uparrow\downarrow$  & $+\dfrac{\hat m_q}{\sqrt{2({\hat m_q}^2 + \textbf{k}_\perp^2)}}$ &
~~~~~~~~~~~~~$\downarrow\uparrow$  & $ - \dfrac{\hat m_q}{\sqrt{2({\hat m_q}^2 + \textbf{k}_\perp^2)}}$ \\
\hline
\end{tabular}
\end{table}

In deriving the complete expression for $D\to a_0(980)$ TFFs, one can take the following correlator
\begin{align}
\Pi_\mu(p,q) = i\int d^4x e^{iq\cdot x}\langle a_0|T \{J_n(x),j_n^\dag(0)\} |0\rangle
\label{ope}
\end{align}
with $J_n(x)=\bar q_1(x)\gamma_\mu\gamma_5 c(x)$, $j_n^\dag(0)=\bar{c}i\gamma_5q_2(0)$. Followed by the standard LCSR approach, one can make the OPE near the light-cone in the space-like region, and insert a complete set of $a_0(980)$-meson states in the physical region. After performing the Borel transformation, we can get the $D\to a_0(980)$ TFFs $f_{\pm}^{D\to a_0}(q^2)$ up to twist-3 accuracy, which can be read off
\begin{align}
  & f_+^{D\to a_0}(q^2)\!=\!\frac{m_c f_{a_0}}{m_D^2 f_D}\int^1_{u_0} du e^{(m_{a_0}^2- s(u))/ M^2} \bigg\{\!-\!\frac{m_c}u \phi_{2;a_0}(u)
\nonumber\\
  &\quad + m_{a_0}\phi_{3,a_0}^p(u)+\frac{m_{a_0}}{6} \bigg[\frac2u \, \phi_{3,a_0}^\sigma(u)-\frac1{m_c^2+u^2m_{a_0}^2-q^2}
\nonumber\\
  &\quad \times \bigg((m_c^2 -u^2m_{a_0}^2+q^2) \, \frac{d\phi_{3,a_0}^\sigma(u)}{du} -\frac{4um_c^2m_{a_0}^2}{m_c^2+ u^2m_{a_0}^2-q^2}
\nonumber\\
  &\quad \times \phi_{3,a_0}^\sigma(u)\bigg)\bigg]\bigg\}
\label{Eq:Da0-TFF}\\
  &f_-^{D\to a_0}(q^2)=\frac{~m_c~ f_{a_0}}{m_D^2 f_D}~\int^1_{u_0} du ~e^{(m_{a_0}^2- s(u))/ M^2}~\bigg[\,\frac{\phi_{3;a_0}^p(u)}{u}
\nonumber\\
  &\quad +\frac{1}{6u}\frac{d\phi_{3,a_0}^\sigma(u)}{du}\bigg].
\label{Eq:Da0-TFF1}
\end{align}
In which the abbreviation $s(u)=(m_b^2+u\bar{u}m_{a_0}^2-\bar{u}q^2)/u$ is used. The lower limits of the integration is $u_0=\{[(s_0 -q^2 -m_{a_0}^2)^2 +4m_{a_0}^2 (m_c^2-q^2)]^{1/2}-(s_0-q^2-m_{a_0}^2)\} / (2m_{a_0}^2)$. Here, $m_D$ and $f_D$ are the mass and decay constant of $D$-meson, $m_c$ is the $c$-quark mass, $\bar u = (1-u)$, $s_0$ stands for the continuum threshold. $\phi_{2;a_0}$ and $\phi_{3,a_0}^\sigma$, $\phi_{3;a_0}^p$ are twist-2 and twist-3 LCDAs, respectively. The twist-3 LCDA can be found in~Refs.~\cite{Han:2013zg,Lu:2006fr}. Then, the decay widths and branching ratio for the semileptonic decay $D\to a_0(980)\ell\bar\nu_\ell$ can be written as~\cite{Cheng:2017fkw,Huang:2021owr}
\begin{align}
  &\frac{d\Gamma(D\to a_0(980)\ell\bar\nu_\ell)}{d q^2} = \frac{G_F^2 |V_{cb}|^2}{768 \pi^3 m_D^3} \frac{(q^2\!-\! m_\ell^2)^2}{q^6}\Big[(m_D^2 \!+\! m_{a_0}^2
\nonumber\\
  &~~-q^2)^2-4m_D^2m_{a_0}^2\Big]^{1/2} \bigg\{\bigg[(q^2+m_{a_0}^2-m_D^2)^2(q^2+2m_{\ell}^2)
\nonumber\\
  &~~-q^2m_{a_0}^2(4q^2 \!+\! 2m_{\ell}^2)\bigg](f_+^{D\to a_0}(q^2))^2 \!+\! 6 q^2 m_{\ell}^2 (m_D^2 -m^2_{a_0}
\nonumber\\
  &~~-q^2) f_+^{D\to a_0}(q^2)f_-^{D\to a_0}(q^2)+6q^4m_{\ell}^2 (f_-^{D\to a_0}(q^2))^2 \bigg\},
\label{wideth}
\nonumber\\
\end{align}
where $m_\ell$, $|V_{cd}|$ and $G_F$ stand for leptonic mass, CKM matrix element and fermi coupling constant.

\section{Numerical analysis}\label{Sec:III}

To do the numerical analysis, we adopt meson's mass $m_{a_0(980)} = 0.980 \pm 0.020~{\rm GeV}$, $m_D=1.865~{\rm GeV}$ and $m_{D^0}=1.870~{\rm GeV}$. Besides, the current quark-mass are $m_u = 2.16^{+0.49}_{-0.26}~{\rm MeV}$ and $m_d = 4.67^{+0.48}_{-0.17}~{\rm MeV}$ at scale $\mu=2~{\rm GeV}$. The $a_0(980)$-meson decay constant is $f_{a_0}=0.409^{+0.022}_{-0.023}~{\rm GeV}$ and mass of the current charm-quark is $m_c=1.237\pm0.02~{\rm GeV}$ at scale $\mu=1.4~{\rm GeV}$. The values of the non-perturbative vacuum condensates appearing in the BFTSR are given below as~\cite{Narison:2014ska, Colangelo:2000dp, Zhong:2021epq},
\begin{eqnarray}
\langle \bar qq\rangle &=&( -2.417_{-0.114}^{+0.227} ) \times 10^{-2}~{\rm GeV}^3
\nonumber\\
\langle g_s\bar q\sigma TGq\rangle &=&(-1.934^{+0.188}_{-0.103}) \times 10^{-2}~{\rm GeV}^5
\nonumber\\
\langle g_s\bar qq\rangle ^2&=&(2.082^{+0.734}_{-0.697}) \times 10^{-3} ~{\rm GeV}^6
\nonumber\\
\langle g_s^2\bar qq\rangle ^2&=&(7.420^{+2.614}_{-2.483}) \times 10^{-3} ~{\rm GeV}^6
\nonumber\\
\langle \alpha_s G^2\rangle&=&0.038 \pm 0.011 ~{\rm GeV}^4
\nonumber\\
\langle g_s^3fG^3\rangle&\simeq& 0.045~{\rm GeV}^6
\nonumber\\
\kappa &=& 0.74\pm0.03
\end{eqnarray}
where the quark-gluon mixed condensate is $\langle g_s\bar q\sigma TGq\rangle = m_0^2 \langle \bar qq\rangle$ with $m_0^2 = 0.80 \pm{0.02}~{\rm GeV}^2$.
The values of the double-quark condensate $q\bar{q}$, quark-gluon mixed condensate $\langle g_s\bar q\sigma TGq\rangle$, four-quark
condensate $\langle g_s^2\bar qq\rangle ^2$, double-gluon condensate $\langle \alpha_s G^2\rangle$ and triple-gluon condensa $\langle g_s^3fG^3\rangle$ are
at $\mu = 2~{\rm GeV}$. Otherwise, these parameters can be calculated to any scale according to the evolution equation.

In the context of BFTSR, there are two important parameters the continuous threshold $s_{a_0}$ and the Borel parameter $M^2$, respectively. We generally take the scale $\mu = M$. Under the 3-loop approximate solution, the $\Lambda_{\rm QCD}^{(n_f)} \simeq (324, 286, 207)~{\rm MeV}$ for the number of quark flavors $n_f = 3, 4, 5$, with which the $ \alpha_s(M_z)=0.1179(10)$, $m_c (\bar{m_c}) = 1.27(2)~{\rm GeV}$, $m_b(\bar{m_b})=4.18^{+0.03}_{-0.02}~{\rm GeV}$ and $M_Z = 91.1876(21)~ {\rm GeV}$ are used. Furthermore, the gluon or quark vacuum condensates and non-perturbative matrix element at initial scale can be running to other scales thought out the renormalization group equations RGE)~\cite{Yang:1993bp, Hwang:1994vp}, which can be written as a general formula:
\begin{align}
\chi(\mu) = \bigg[\frac{\alpha_s(\mu_0)}{\alpha_s(\mu)}\bigg]^{y(n_f)} \chi(\mu_0)
\end{align}
with which the function $y(n_f)$ are $-4/b$, $-4/b$ and $-2/(3b)$ for the $m_q$, $\langle \bar qq\rangle$ and $\langle g_s\bar q\sigma TGq\rangle$, respectively. The coefficient $b = (33-2n_f)/3$ and $n_f$ is the number of active quark flavors. According to the basic assumption of BFTSR, it is worth noting that $g_s$ is the coupling constant between background fields in the above vacuum condensates, which is different from the coupling constant in pQCD and should be absorbed into the vacuum condensates as part of these non-perturbative parameters. The RGE of the Gegenbauer moments of the $a_0(980)$-meson twist-2 LCDA is ~\cite{Cheng:2005nb}:
\begin{eqnarray}
a_n^{2;a_0}(\mu) &=& a_n^{2;a_0}(\mu_0) E_n(\mu, \mu_0),
\label{Eq:anRGE}
\end{eqnarray}
with $E_n(\mu, \mu_0) = [\alpha_s(\mu)/\alpha_s(\mu_0)]^{-(\gamma_n^{(0)}+4)/(b)}$. The $\mu_0$ and $\mu$ represent the initial scale and the running scale, the one-loop anomalous dimensions is $\gamma_n^{(0)} = C_F \{[1-2/[(n+1)(n+2)]+4\sum^{n+1}_{j=2}1/j\}$, with $C_F = 4/3$. According to the RGE of the Gegenbauer moments, one can get the moments $\langle\xi^n_{2;a_0}\rangle|_\mu$ for the arbitrary scale $\mu$.
\begin{table}[t]
\centering
\renewcommand\arraystretch{1.3}
\small
\caption{The determined Borel windows and the corresponding $a_0(980)$-meson twist-2 LCDA moments $\langle\xi^n_{2;a_0}\rangle|_\mu$ with $n=(1,3,5,7,9)$ at the scale $\mu_k= 1.4~{\rm GeV}$. Where all input parameters are set to be their central values. In which the abbreviation ``Con.'' indicate the continuum contributions.}
\begin{tabular}{ l c c l}
\\ \hline
$n$~~~~~~& ~~~~~~~~$M^2$~~~~~~~~&~~~~~~~~~~~~$\langle\xi^n_{2;a_0}\rangle|_\mu$~~~~~~~~~~~~&~~~Con. \\
\hline
$1$~ & ~$[2.483,3.483]$~ & ~$[-0.250,-0.203]$~& ~$<30\%$~\\
$3$~ & ~$[1.869,2.869]$~ & ~$[-0.095,-0.129]$~& ~$<20\%$~\\
$5$~ & ~$[3.057,4.057]$~ & ~$[-0.041,-0.074]$~& ~$<40\%$~\\
$7$~ & ~$[4.143,5.143]$~ & ~$[-0.020,-0.048]$~& ~$<60\%$~\\
$9$~ & ~$[4.916,5.916]$~ & ~$[-0.000,-0.027]$~& ~$<75\%$~\\
\hline
\end{tabular}
\label{Tab:tbw}
\end{table}

Furthermore, the continuum threshold parameter for the sum rule $\langle\xi^n_{2;a_0}\rangle|_\mu$ can be determined by normalization for $\langle\xi^{p;0}_{3;a_0}\rangle|_\mu$, which leads to $s_{a_0}=7~{\rm GeV}^2$. Then, the Borel window for the each order of $a_0(980)$-meson LCDA moments can be determined by limiting the continuum states and the dimension-six condensates contributions. Then, the moments $\langle\xi^n_{2;a_0}\rangle|_\mu$ with $n=(1,3,5,7,9)$ within uncertainties coming from Borel parameters are listed in Table~\ref{Tab:tbw}. Based on the BFTSR, the dimension-six condensates contributions for $\langle\xi^n_{2;a_0}\rangle|_\mu$ are less than $1\%$ for all the $n$th-order. To get the suitable Borel window, the continuum contributions for $\langle\xi^n_{2;a_0}\rangle|_\mu$ are restrict to $(30, 20, 40, 60, 75)\%$ for $n=(1,3,5,7,9)$ respectively. Since the dimension-six condensates contributions are very small, we determine the upper limit of the Borel parameter $M^2$ through the continuum contributions. Then the lower limit of the Borel parameter $M^2$ can be determined by the method of the upper limits, so as to obtain the appropriate Borel window. At the same time, the values of the moments $\langle\xi^n_{2;a_0}\rangle|_\mu$ are stable in the appropriate Borel window. To have a deeper insight into the relationship of the LCDA moments versus Borel parameter $M^2$, the first five moments's curves are shown in Fig.~\ref{Fig:fxinM2}, which can be seen that
\begin{itemize}
 \item In the region $M^2 \in {\rm{[1,2]}}$ of Borel window, the curves for $\langle\xi^n_{2;a_0}\rangle|_\mu$ changed dramatically. With the increase of Borel window $M^2$, the change trend of moments $\langle\xi^1_{2;a_0}\rangle|_\mu$, $\langle\xi^3_{2;a_0}\rangle|_\mu$,
   $\langle\xi^5_{2;a_0}\rangle|_\mu$,
   $\langle\xi^7_{2;a_0}\rangle|_\mu$, $\langle\xi^9_{2;a_0}\rangle|_\mu$ tends to be gentle.
 \item With $n$ increases, the absolute value of the moments $\langle\xi^n_{2;a_0}\rangle|_\mu$ tend to be smaller.
\end{itemize}
\begin{figure}[t]
\centering
\includegraphics[width=0.45\textwidth]{Fig1_Borel.eps}
\caption{(Color online) The $a_0(980)$-meson leading-twist DA moments $\langle\xi^n_{2;a_0}\rangle|_\mu$ with $n=(1,3,5,7,9)$ versus the Borel parameter $M^2$, where all input parameters are set to be their central values. In which, the shaded band indicate the Borel Windows for $n=(1,3,5,7,9)$, respectively.}
\label{Fig:fxinM2}
\end{figure}
\begin{figure*}[t]
\centering
\includegraphics[width=0.45\textwidth]{Fig2a_DA.eps}
\includegraphics[width=0.45\textwidth]{Fig2b_DA_mq.eps}
\caption{(Color online) The $a_0(980)$-meson twist-2 DA curves in this work. As a comparison, we present the curve trend of QCDSR distribution amplitude~\cite{Cheng:2005nb} in the left panel. In the right panel, the $a_0(980)$-meson twist-2 DA $\phi_{2;a_0}(x,\mu)$ at $\mu_k = 1.4~{\rm GeV}$ with the constituent quark mass $m_q = (200, 250, 300, 350)~{\rm MeV}$ are given, respectively.}
\label{fDAcures}
\end{figure*}
Taking all the input uncertainties into consideration, we can get the moments $\langle\xi^n_{2;a_0}\rangle|_\mu$ with $n=(1,3,5,7,9)$ under two types of factorization scale $\mu_0$ and $\mu_k$,
\begin{align}
&\langle\xi^1_{2;a_0}\rangle|_{\mu_0} = -0.307(43), && \langle\xi^1_{2;a_0}\rangle|_{\mu_k} = -0.247(34), \nonumber\\
&\langle\xi^3_{2;a_0}\rangle|_{\mu_0} = -0.181(34), && \langle\xi^3_{2;a_0}\rangle|_{\mu_k} = -0.123(22), \nonumber\\
&\langle\xi^5_{2;a_0}\rangle|_{\mu_0} = -0.078(28), && \langle\xi^5_{2;a_0}\rangle|_{\mu_k} = -0.065(21), \nonumber\\
&\langle\xi^7_{2;a_0}\rangle|_{\mu_0} = -0.049(26), && \langle\xi^7_{2;a_0}\rangle|_{\mu_k} = -0.040(19), \nonumber\\
&\langle\xi^9_{2;a_0}\rangle|_{\mu_0} = -0.036(24), && \langle\xi^9_{2;a_0}\rangle|_{\mu_k} = -0.019(16),
\label{Eq:xin_mu0}
\end{align}
Our results for the first two order are slightly smaller than the Cheng's predictions, e.g. $\langle\xi^1_{2;a_0}\rangle|_{\mu_0} = -0.56(5)$ and $\langle\xi^3_{2;a_0}\rangle|_{\mu_0} = -0.21(3)$ by using the QCDSR approach from Ref.~\cite{Cheng:2005nb}, which is more likely to be antisymmetric behavior. The little difference may be related to the different methods in determining the continuum threshold $s_{a_0}$. Meanwhile, the higher order such as $n=(5,7,9)$ are given for the first time.

Secondly, the free parameters $\alpha_{2;a_0}$ and $A_{2;a_0}$ for the LCHO model can be fixed by adopting the method of least squares to fit $\langle\xi^n_{2;a_0}\rangle|_\mu$ calculated in the framework of the BFTSR shown in Eq.~\eqref{Eq:xin_mu0}. The goodness of fit can be judged by the probability $P_{\chi^2}$ ($P_{\chi^2} \in [0,1]$)~\footnote{For more details, one can see our previous work for pion~\cite{Zhong:2021epq}.}.
By fitting the moments $\langle\xi^n_{2;a_0}\rangle|_{\mu_k}$ at the scale $\mu_k=1.4~{\rm GeV}$, the optimal model parameters are obtained as follows,
\begin{eqnarray}
A_{2;a_0}     &=& -300 {\rm GeV}^{-1}, \nonumber\\
\alpha_{2;a_0}  &=& -0.25,         \nonumber\\
\beta_{2;a_0}   &=& 0.5 {\rm GeV},
\label{ModelParameterIII}
\end{eqnarray}
The above parameters are derived from $m_q=250~{\rm MeV}$. From these parameters, we obtain the scalar $a_0(980)$-meson twist-2 LCDA $\phi_{2;a_0}(x,\mu_k)$, which is shown in the left panel of Fig.~\ref{fDAcures}. In the meantime, we also present the results of QCDSR~\cite{Cheng:2005nb} prediction as a comparison. From the perspective of variation trend, our prediction results are relatively consistent with the those of QCDSR, and twist-2 LCDA are both antisymmetric. However, there are also some differences between the two, which may be due to differences in calculation methods. The former uses the Gegenbauer moments polynomials expansion, while we take the LCHO model to construct LCDA $\phi_{2;a_0}(x,\mu)$.
\begin{table}
\centering
\renewcommand\arraystretch{1.3}
\small
\caption{The model parameters $m_q$ (in unit: MeV) of LCHO model $\varphi_{2;a_0}^{\rm IV}(x)$ under different quark masses and their corresponding goodness of fit.}\label{table:fitting}
\begin{tabular}{lllll}
\\ \hline
$m_q$~~~~~~~~~~~~~~~~~~ & $A_{2;a_0}(\rm GeV^{-1})$ & $\alpha_{2;a_0}$~~~~~~~~~~~~~~ & $P_{\chi^2_{\rm min}}$ \\ \hline
200MeV & $-440$ & $-0.08$ & $0.912$\\
250MeV & $-300$ & $-0.25$ & $0.920$\\
300MeV & $-190$ & $-0.48$ & $0.928$\\
350MeV & $-140$ & $-0.66$ & $0.935$ \\
 \hline
\end{tabular}
\end{table}
\begin{figure}[t]
\centering
\includegraphics[width=0.45\textwidth]{Fig3_Fitting.eps}
\caption{(Color online) Relationship between the goodness of fit $P_{\chi^2_{\rm min}}$ and the two LCHO parameters $A_{2;a_0}$, $\alpha_{2;a_0}$. In which the $P_{\chi^2_{\rm min}}$ is separated into three area with different color.}
\label{Fig:curves}
\end{figure}
\begin{figure*}
\centering
\includegraphics[width=0.45\textwidth]{Fig4a_fp.eps}\includegraphics[width=0.45\textwidth]{Fig4b_fm.eps}
\caption{(Color online) The $D\to a_0(980)\ell\bar{\nu}_{\ell}$ TFFs $f_+^{D\to a_0}(q^2)$ and $f_-^{D\to a_0}(q^2)$ of our predictions within uncertainties. In which, the darker band stand for the results calculated in LCSR, and the light colored band stand for the SSE. Meanwhile, the results for LCSR~\cite{Huang:2021owr} and Ads/QCD~\cite{Momeni:2022gqb} as a comparison.}
\label{Fig:TFF}
\end{figure*}

For further study, we also analyze the relationship between the goodness-fit of the $a_0 (980)$-meson twist-2 LCDA $\phi_{2;a_0}(x,\mu_k)$ and the quark mass $m_q$ listed in Table~\ref{table:fitting}. It is obvious that as the quark mass increases, the goodness of fit also increases. The goodness of fit $P_{\chi^2_{\rm min}}$ also indicates that $\beta_{2;a_0}=0.5 {\rm GeV}$ is reasonable. Considering the effect of quark mass $m_q$ on the $a_0(980)$-meson twist-2 LCDA $\phi_{2;a_0}(x,\mu_k)$, the behavior of LCHO model $\phi_{2;a_0}(x,\mu)$ for different quark masses $m_q=(200,250,300,350)~{\rm MeV}$ is shown in the right panel of Fig.~\ref{fDAcures}. As can be seen from the figure, the peak value of the LCHO model curves increases with the increase of mass $m_q$. In addition to this, the relationship between goodness of fit $P_{\chi^2_{\rm min}}$ and parameters $A_{2;a_0}$ and $\alpha_{2;a_0}$ is also shown in Fig.~\ref{Fig:curves}. From the table we can see that the $P_{\chi^2_{\rm min}}$ can reach to $93\%$, which shows the fitting is good.

The TFFs is an important parameter in the calculation of the semileptonic decay $D\to a_0(980)\ell\bar \nu_\ell$. To obtain the numerical results of TFFs, we take $\mu_k = 1.4~{\rm GeV}$, $m_D = 1.869~{\rm GeV}$~\cite{PDGnew}, $f_{a_0(980)} = (0.409^{+0.022}_{-0.023})~{\rm GeV}$~\cite{{Cheng:2005nb}}. For the continuum threshold $s_0$, it is generally taken near the mass square of the first excited state of $D$-meson, that is, near the mass square of $D(2550)$-meson. Based on the prediction of the sum rule of heavy quark effective theory (HQET)~\cite{Huang:1998sa}, we take the continuum threshold parameter $s_0= 6.50(25)~{\rm GeV}^2$, which is also agree with Ref.~\cite{Li:2012gr}.
\begin{table}[t]
\centering
\renewcommand\arraystretch{1.3}
\caption{TFFs at the point of large recoil $q^2\simeq 0$ for $D\to a_0(980)$ within uncertainties. And a comparison with other theoretical groups are also given.}\label{factor}.
\small
\begin{tabular}{l l l}
\hline
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~&$f_+^{D\to a_0}(0)$~~~~~~~~~~~~~~~~~~& $f_-^{D\to a_0}(0)$ \\
\hline
This~work & $1.058^{+0.068}_{-0.035}$ &	$0.764^{+0.044}_{-0.036}$  \\
CCQM~\cite{Soni:2020sgn} & $0.55^{+0.02}_{-0.02}$ & $0.03^{+0.01}_{-0.01}$ \\
LCSR 2021~\cite{Huang:2021owr} &$0.85^{+0.10}_{-0.11}$	& $\!\!\!\!\!-0.85^{+0.10}_{-0.11}$ \\
LCSR 2017~\cite{Cheng:2017fkw} & $1.76(26)$ & $0.31(13)$ \\
 AdS/QCD~\cite{Momeni:2022gqb}&$0.72(9)$&~~~~-\\
 \hline
\end{tabular}
\end{table}

According to LCSR, in the suitable Borel window, we predict the end values of the $D\to a_0(980)\ell\bar{\nu}_{\ell}$ semileptonic decay TFFs shown in Table~\ref{factor}. As a comparison, the results predicted from various approaches, CCQM ~\cite{Soni:2020sgn}, LCSR~2017~\cite{Cheng:2017fkw}, LCSR~2021~\cite{Huang:2021owr} and AdS/QCD~\cite{Momeni:2022gqb} also are presented in Table~\ref{factor}. It is not difficult to see from the Table~\ref{factor} that our predicted results are significantly different from those obtained by other predictions. The reason lies in the twist-2 LCDA $\phi_{2;a_0}(x,\mu)$ is different, which the LCHO model is used here.

\begin{figure*}[t]
\centering
\includegraphics[width=0.45\textwidth]{Fig51_D0e.eps}~~\includegraphics[width=0.45\textwidth]{Fig52_D0mu.eps}
\includegraphics[width=0.45\textwidth]{Fig53_Dpe.eps}~~\includegraphics[width=0.45\textwidth]{Fig54_Dpmu.eps}
\caption{(Color online) The differential decay widths for $D\to a_0(980)\ell\bar{\nu}_{\ell}$ with $\ell = (e,\mu)$ decay process within uncertainties. In which, the LCSR~\cite{Huang:2021owr} and AdS/QCD~\cite{Momeni:2022gqb} predictions are also presented.}
\label{width}
\end{figure*}
In the physical sense, the LCSR method is suitable for in low and intermediate $q^2$ region. The physical region for $D\to a_0(980)$ TFFs is $m_e\le q^2\le 0.48~{\rm GeV}^2$. In order to obtain reasonable LCSR results, we can extend them to the entire physical $q^2$-region $m_e\le q^2\le (m_D^2-m_{a_0})^2 = 0.8~{\rm GeV}^2$. Therefore, we can fit the complete analysis results with simplified series expansion (SSE), which is a rapidly convergent series on the $z(t)$-expansion~\cite{Fu:2018yin,Bharucha:2015bzk,Bourrely:2008za}.
\begin{align}
{{f}_{i}}(q^2)={P_i}(q^2)\sum\limits_{{k = 0,1,2}} a_k^i {{[z(}{q^2}) - z(0)]}^{k},
\label{ff}
\end{align}
where ${{f}_{i}}(q^2)$ stands for $D\to a_0(980)\ell\bar{\nu}_{\ell}$ TFFs, $a_k^i$ the fit coefficients,and
\begin{align}
P_i(q^2)=\bigg(1-\frac{q^2}{m_{R,i}^2}\bigg)^{-1}
\nonumber\\
z(t)=\frac{\sqrt{t_+-t}-\sqrt{t_+-t_0}}{\sqrt{t_+-t}+\sqrt{t_+-t_0}}
\end{align}
\begin{table}[t]
\centering
\renewcommand\arraystretch{1.3}
\caption{The fitting parameters $a_i$ with $i=(1,2)$ for TFFs $f^{D\to a_0}_+(q^2)$ and $f^{D\to a_0}_-(q^2)$. In which, the goodness of fit $\Delta$ are also present.}\label{Tab:SSE}
\small
\begin{tabular}{l l l l}
\hline
$f^{D\to a_0}_+(q^2)$~~~   & Central value~~~ & Upper limits~~~ & Lower limits \\ \hline
$a_1$     & $7.279$  & $6.686$  & $7.562$  \\
$a_2$     & $33.88$ & $38.93$ & $30.14$  \\
$\Delta$    & 0.19\textperthousand & 0.19\textperthousand & 0.19\textperthousand \\ \hline
$f^{D\to a_0}_-(q^2)$~~~   & Central value~~~ & Upper limits~~~ & Lower limits \\ \hline
$a_1$    & \!\!\!\!\!$-1.556$  & \!\!\!\!\!$-1.776$  & \!\!\!\!\!$-1.178$ \\
$a_2$    & $147.1$ & $167.9$ & $129.1$ \\
$\Delta$  & 0.29\textperthousand & 0.33\textperthousand & 0.25\textperthousand \\
\hline
\end{tabular}
\end{table}
with $t_{\pm}\equiv(m_D\pm m_{a_0})^2$, $t_0=t_+(1-\sqrt{1- t_-/t_+})$. $P_i(q^2)$ is a simple pole corresponding to the first resonance in the spectrum. The fitting parameters $a_i$ for the TFFs within uncertainties are given in Table~\ref{Tab:SSE}. Meanwhile, the goodness of fit $\Delta = {\sum_t|F_i(t) - F_i^{\rm fit}(t)|}/{\sum_t|F_i(t)|} \times 100$ with $t\in[0,1/100,\cdots,100/100]\times 0.48~{\rm GeV}$ for each TFFs are also presented. After extrapolating TFFs to the whole $q^2$ region, the curves $f_{\pm}^{D\to a_0}(q^2)$ are shown in the Fig.~\ref{Fig:TFF}. We also give the curves obtained by LCSR 2021~\cite{Huang:2021owr} and AdS/QCD~\cite{Momeni:2022gqb} for comparison. The results show our prediction is a certain discrepancy between the results in Ref.~\cite{Huang:2021owr}. There are also some deviations from results predicted by AdS/QCD, but the curve trend is relatively consistent. In most cases, The $f_+^{D\to a_0}(q^2)$ shows an upward trend with the increase of $q^2$. Our prediction is quite reasonable. The $f_-^{D\to a_0}(q^2)$ is also exhibit upward tendency and it turns out that there are some differences.

As other important parameters for $D\to a_0(980)\ell\bar{\nu}_{\ell}$, the CKM matrix element $|V_{cd}|= 0.221\pm 0.008$ from the PDG~\cite{PDGnew}, the fermi coupling constant $G_F=1.166\times10^{-5}~\rm{GeV}^{-2}$. By taking the derived $D\to a_0(980)\ell\bar{\nu}_{\ell}$ TFFs and the related parameters into the differential decay widths Eq.~\eqref{wideth}, one can get the differential decay widths of $D\to a_0(980)\ell\bar{\nu}_{\ell}$ with $\ell= (e,\mu)$ presented in Fig.~\ref{width}. For comparison, we also present the results of LCSR 2021~\cite{Huang:2021owr} predictions in Fig.~\ref{width}. The figure shows that our prediction is consistent with the result of LCSR 2021 within a certain error range. In Fig.~\ref{width}(a), discrepancy between our prediction and the result of AdS/QCD, which may be caused by the different TFFs. Obviously, the predicted results converge to zero in the small recoil region $q^2=(m_D-m_{a_0})^2$. The shaded part in the figure is the error of the width, which mainly comes from the uncertainty of all parameters.

\begin{table*}[t]
\centering
\renewcommand\arraystretch{1.3}
\small
\caption{Branching fractions for the four different channels of $D\to a_0(980)\ell\bar{\nu}_{\ell}$ (in unit: $10^{-4}$). To make a comparison, we also listed the CCQM~\cite{Soni:2020sgn}, LCSR~\cite{Huang:2021owr,Cheng:2017fkw} and AdS/QCD~\cite{Momeni:2022gqb} predictions.}
\begin{tabular}{lllll}
\hline
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~&$ D^0\to {a_0(980)}^{-}e^+\nu_e$~~~~~~~~&$D^0\to {a_0(980)}^{-}{\mu}^+\nu _\mu$~~~~~~~~&$D^+\to {a_0(980)}^0e^+\nu_e$~~~~~~~~&$D^+\to {a_0(980)}^0{\mu}^+\nu_\mu$ \\ \hline
This work    &$1.529^{+0.230}_{-0.108}$ & $1.452^{+0.220}_{-0.105}$& $1.927^{+0.288}_{-0.137}$&$1.830^{+0.277}_{-0.130}$ \\
CCQM ~\cite{Soni:2020sgn} & $1.68\pm{0.15}$  &$1.63\pm{0.14}$  & $2.18\pm{0.38}$ & $2.12\pm{0.37}$ \\
LCSR 2017~\cite{Cheng:2017fkw} &$4.08^{+1.37}_{-1.22}$   & -	& $5.40^{+1.78}_{-1.59}$ &-\\
LCSR 2021~\cite{Huang:2021owr}  &$1.36$   & $1.21$	& $1.79$ &$1.59$\\
AdS/QCD~\cite{Momeni:2022gqb}&$2.44\pm0.30$&-&-&-\\
\hline
\label{branching ratios}
\end{tabular}
\end{table*}
\begin{table*}[t]
\centering
\renewcommand\arraystretch{1.3}
\small
\caption{The absolute branching ratio of $D\to a_0(980)(\to\eta\pi)e^+\nu_e$ (unit: $10^{-4}$) of our predictions within uncertainties. To make a comparison, we also listed the BESIII collaboration~\cite{BESIII:2018sjg}, LCSR results~\cite{Huang:2021owr} and PDG average value~\cite{PDGnew}.}
\begin{tabular}{lllll}
\hline
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~&${\cal B}(D^0\to a_0(980)^-(\to\eta \pi^-) e^+\nu_e)$~~~~~~~~~&${\cal B}(D^+ \to a_0(980)^0(\to\eta \pi^0) e^+\nu_e)$ \\ \hline
This work     & $1.292^{+0.194}_{-0.092}$& $1.628^{+0.244}_{-0.116}$ \\
BESIII~\cite{BESIII:2018sjg}        & $1.33^{+0.33}_{-0.29}$ & $1.66^{+0.81}_{-0.66}$ \\
LCSR 2021~\cite{Huang:2021owr}       & $1.15$	& $1.51$ \\
PDG~\cite{PDGnew}        & $1.33^{+0.30}_{-0.29}$ & $1.7^{+0.8}_{-0.7}$ \\
\hline
\label{ratios}
\end{tabular}
\end{table*}
Furthermore, by taking the lifetimes of $D^0$, $D^+$ meson $\tau_{D^0}=(0.410\pm0.001) \rm {ps}$ and $\tau _{D^+}=(1.033\pm0.005) \rm {ps}$ from PDG, we can obtain the branching ratio of semileptonic decay channels $D^+\to {a_0(980)}^0\ell^+\nu_\ell$ ($\ell= e,\mu$) and $D^0\to {a_0(980)}^-\ell^+\nu_\ell$ in Table~\ref{branching ratios}. What's more, the predictions from theoretical groups such as CCQM~\cite{Soni:2020sgn}, LCSR~\cite{Cheng:2017fkw,Huang:2021owr}, AdS/QCD~\cite{Momeni:2022gqb} are shown in Table \ref{branching ratios}. It is obvious that our results are relatively consistent with those of CCQM and LCSR 2021 within the error. However, the significant difference between our predictions and LCSR 2017, AdS/QCD may be caused by the difference of the TFFs and $a_0(980)$-meson twist-2 LCDA.

Additionally, we also calculate the absolute branching ratio of decays $D\to \eta\pi\ell\bar{\nu}_{\ell}$ by using the relationship
\begin{align}
&{\cal B}(D\to {a_0(980)} (\to \eta\pi)e^+\nu_e) = \nonumber\\
&\quad{\cal B}(D\to {a_0(980)}e^+\nu_e)\times{\cal B}(a_0(980)\to \eta \pi).
\end{align}
In which the ${\cal B}(a_0(980)^0 \to \eta \pi^0)={\cal B}(a_0(980)^- \to \eta \pi^-) = 0.845 \pm0.017$ can be used~\cite{Cheng:2005nb}. Combing with the ${\cal B}(D\to {a_0(980)}e^+\nu_e)$ been calculated in this paper, we can get the results for ${\mathcal B}(D^0\to a_0(980)^-(\to \eta \pi^-) e^+\nu_e)$ and ${\mathcal B}(D^+\to a_0(980)^0(\to \eta \pi^0) e^+\nu_e)$, which are listed in Table~\ref{ratios}. For comparison, the results of the BESIII~\cite{BESIII:2018sjg} collaboration, theory groups LCSR 2021~\cite{Huang:2021owr} and PDG~\cite{PDGnew} predictions are also given. The results show that our predicted absolute branching ratio are in good agreement with those predicted of BESIII, LCSR 2021 and PDG within the error range. It is clear that our prediction is more accurate than the LCSR 2021 prediction, with an improvement of about $12\%$. It shows that the result of our prediction is reasonable.

\section{Summary}\label{Sec:IV}

In this paper, we have calculated the moments of $a_0(980)$-meson twist-2 LCDA by adopting the QCDSR approach within the background field theory. The continuum threshold parameter $s_{a_0}$ is determined from the normalization for the $a_0(980)$-meson twist-3 LCDA $0_{\rm th}$-order moments. After seeking the suitable Borel windows, we present the first five order moments, i.e. $\langle\xi^n_{2;a_0}\rangle|_\mu$ with $n = (1,3,5,7,9)$ under two different factorization scale $\mu_0$ and $\mu_k$. Then, we study $a_0(980)$-meson twist-2 LCDA $\phi_{2;a_0}(x,\mu)$ based on the LCHO model for improving the accuracy of the calculation. The least square method is used to fit the moments $\langle\xi^n_{2;a_0}\rangle|_\mu$ and to determine the model parameters. The goodness of fit can be up to $93\%$. Then, the curves of $a_0(980)$-meson twist-2 LCDA comparing with other theoretical groups and with different constituent quark masses are presented.

The $D\to a_0(980)$ TFFs are calculated within the LCSR approach. The TFFs at large recoil region are listed in Table~\ref{factor}. After extrapolating the TFFs to the whole $q^2$ region, the curves $f_{\pm}^{D\to a_0}(q^2)$ are shown in Fig.~\ref{Fig:TFF}. A comparison of TFFs with other LCSR and AdS/QCD predictions are also given. Using the resultant TFFs, we further studied the semileptonic decays $D\to a_0(980)\ell\bar{\nu}_{\ell}$ with $\ell = (e,\mu)$. Their differential decay widths are presented in Fig.~\ref{width}, and their branching fractions are given in Table~\ref{branching ratios}. The ratio of partial branching fractions is given
\begin{align}
\dfrac{{\cal B}(D^0\to a_0(980)^- e^+ \nu_e)}{{\cal B}(D^+\to a_0(980)^0 e^+ \nu_e)} = 0.793^{+0.180}_{-0.159},
\end{align}
which agree with the CCQM prediction~\cite{Soni:2020sgn} and the LCSR predictions~\cite{Cheng:2017fkw, Huang:2021owr} within errors.

After considering the decay $a_0(980)\to \eta \pi$, we have calculated the branching fractions for the decay processes $D\to \eta\pi e^+\nu_e$, ${\mathcal B}(D^0\to \eta \pi^- e^+ \nu_e) = (1.292^{+0.194}_{-0.092})\times10^{-4}$ and ${\mathcal B}(D^+\to \eta \pi^0 e^+\nu_e)=(1.628^{+0.244}_{-0.116})\times10^{-4}$. These results are consistent with the BESIII data and the PDG average value within errors. Along with other results of branching fraction for scalar meson discovered experimentally, we will have a reliable input for understanding the nature of the light scalar mesons.

\acknowledgments

Hai-Bing Fu and Tao Zhong would like to thank the Institute of Theoretical Physics in Chongqing University for kind hospitality. This work was supported in part by the National Natural Science Foundation of China under Grant No.12265010, No.12265009, No.12175025 and No.12147102, the Project of Guizhou Provincial Department of Science and Technology under Grant No.ZK[2021]024, the Project of Guizhou Provincial Department of Education under Grant No.KY[2021]030, and by the Chongqing Graduate Research and Innovation Foundation under Grant No. ydstd1912.