#  $\eta^{(\prime)}$ -meson twist-2 distribution amplitude within QCD sum rule
approach and its application to the semi-leptonic decay
$D_{s}^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$

Dan-Dan Hu1 Hai-Bing Fu1,3<EMAIL_ADDRESS>(Corresponding author) Tao Zhong1
Long Zeng2,3 Wei Cheng4 Xing-Gang Wu2,3<EMAIL_ADDRESS>1Department of
Physics, Guizhou Minzu University, Guiyang 550025, P.R. China  2Department of
Physics, Chongqing University, Chongqing 401331, P.R. China  3Chongqing Key
Laboratory for Strongly Coupled Physics, Chongqing 401331, P.R. China  4School
of Science, Chongqing University of Posts and Telecommunications, Chongqing
400065, P.R. China

###### Abstract

In this paper, we make a detailed discussion on the $\eta$ and
$\eta^{\prime}$-meson leading-twist light-cone distribution amplitude
$\phi_{2;\eta^{(\prime)}}(u,\mu)$ by using the QCD sum rules approach under
the background field theory. Taking both the non-perturbative condensates up
to dimension-six and the next-to-leading order (NLO) QCD corrections to the
perturbative part, its first three moments
$\langle\xi^{n}_{2;\eta^{(\prime)}}\rangle|_{\mu_{0}}$ with $n=(2,4,6)$ can be
determined, where the initial scale $\mu_{0}$ is set as the usual choice of
$1$ GeV. Numerically, we obtain
$\langle\xi_{2;\eta}^{2}\rangle|_{\mu_{0}}=0.231_{-0.013}^{+0.010}$,
$\langle\xi_{2;\eta}^{4}\rangle|_{\mu_{0}}=0.109_{-0.007}^{+0.007}$, and
$\langle\xi_{2;\eta}^{6}\rangle|_{\mu_{0}}=0.066_{-0.006}^{+0.006}$ for
$\eta$-meson,
$\langle\xi_{2;\eta^{\prime}}^{2}\rangle|_{\mu_{0}}=0.211_{-0.017}^{+0.015}$,
$\langle\xi_{2;\eta^{\prime}}^{4}\rangle|_{\mu_{0}}=0.093_{-0.009}^{+0.009}$,
and
$\langle\xi_{2;\eta^{\prime}}^{6}\rangle|_{\mu_{0}}=0.054_{-0.008}^{+0.008}$
for $\eta^{\prime}$-meson. Next, we calculate the $D_{s}\to\eta^{(\prime)}$
transition form factors (TFFs) $f^{\eta^{(\prime)}}_{+}(q^{2})$ within QCD
light-cone sum rules approach up to NLO level. The values at large recoil
region are $f^{\eta}_{+}(0)=0.476_{-0.036}^{+0.040}$ and
$f^{\eta^{\prime}}_{+}(0)=0.544_{-0.042}^{+0.046}$. After extrapolating TFFs
to the allowable physical regions within the series expansion, we obtain the
branching fractions of the semi-leptonic decay, i.e.
$D_{s}^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$, i.e. ${\cal
B}(D_{s}^{+}\to\eta^{(\prime)}e^{+}\nu_{e})=2.346_{-0.331}^{+0.418}(0.792_{-0.118}^{+0.141})\times
10^{-2}$ and ${\cal
B}(D_{s}^{+}\to\eta^{(\prime)}\mu^{+}\nu_{\mu})=2.320_{-0.327}^{+0.413}(0.773_{-0.115}^{+0.138})\times
10^{-2}$ for $\ell=(e,\mu)$ channels respectively. And in addition to that,
the mixing angle for $\eta-\eta^{\prime}$ with $\varphi$ and ratio for the
different decay channels ${\cal R}_{\eta^{\prime}/\eta}^{\ell}$ are given,
which show good agreement with the recent BESIII measurements.

##  I Introduction

The $D^{+}_{s}$-meson, which is composed of a charm quark and a strange
antiquark, has been discovered in year 1993. There are rich physics contents
in $D^{+}_{s}$-meson decays. The $D_{s}^{+}$-meson semileptonic or exclusive
decay processes provide important heavy-to-light theoretical bases for
studying heavy quark decays, investigating light meson spectroscopy and
supplying a bridge between weak and strong interaction couplings of quarks.
More and more experimental results have been reported from the BABAR, the
CLEO, the BESIII collaborations, and etc., such as the
$D_{s}^{+}\to(\eta,\eta^{\prime},K^{0},a_{0}(980),f_{0}(980),\phi,K^{*0})\ell^{+}\nu_{\ell}$
decays’ branching fractions are within the range of $[0.12,2.61]\%$
BaBar:2008gpr ; CLEO:2009dyb ; CLEO:2009ugx ; BESIII:2017ikf ; BESIII:2018xre
; BESIII:2021duu . Total semileptonic branching fractions provide useful
discrimination on the different theoretical evaluations of hadronic matrix
elements, which sizably affect the charm quark semileptonic decays. The
$D_{s}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$ is different from the usually
considered channels with final state containing light-quark composition only,
and it has attracted much attention from both theoretical and experimental
groups. Moreover, the $\eta^{(\prime)}$-mesons, composed by $s\bar{s}$ quark
pair, are especially intriguing, since the $s$-quark plays an important role
for the flavor physics. In deep leaning of those two processes, one can obtain
useful information on the CKM matrix element $|V_{cs}|$ and the heavy-to-light
transition form factors (TFFs).

Experimentally, the semileptonic decay processes for
$D_{s}^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$ have been found by the CLEO
collaboration early in year 1995, and their measured value of the ratio of
branching fractions ${\cal B}(D_{s}^{+}\to\eta^{\prime}e^{+}\nu_{e})/{\cal
B}(D_{s}^{+}\to\eta e^{+}\nu_{e})$ is $0.35\pm 0.09\pm 0.07$
Brandenburg:1995qq . Then, the CLEO collaboration issued the measured value of
the branching fractions in years 2009 and 2015, i.e. ${\cal
B}(D_{s}^{+}\to\eta e^{+}\nu_{e})=(2.48\pm 0.29\pm 0.13)\%$, ${\cal
B}(D_{s}^{+}\to\eta e^{+}\nu_{e})=(2.28\pm 0.14\pm 0.20)\%$, ${\cal
B}(D_{s}^{+}\to\eta^{\prime}e^{+}\nu_{e})=(0.91\pm 0.33\pm 0.05)\%$, and
${\cal B}(D_{s}^{+}\to\eta^{\prime}e^{+}\nu_{e})=(0.68\pm 0.15\pm 0.06)\%$
Yelton:2009aa ; Hietala:2015jqa . In year 2017, the BESIII collaboration
measured the branching fractions by using the same channels based on the
integrated luminosity of $482~{}{\rm pb^{-1}}$ of the $e^{+}e^{-}$ collision
at the center-of-mass energy $\sqrt{s}=4.009~{}{\rm GeV}$, and they issued
${\cal B}(D_{s}^{+}\to\eta\mu^{+}\nu_{\mu})=(2.42\pm 0.46\pm 0.11)\%$, ${\cal
B}(D_{s}^{+}\to\eta^{\prime}\mu^{+}\nu_{\mu})=(1.06\pm 0.54\pm 0.07)\%$,
${\cal B}(D_{s}^{+}\to\eta e^{+}\nu_{e})=(2.30\pm 0.31\pm 0.08)\%$, and ${\cal
B}(D_{s}^{+}\to\eta^{\prime}e^{+}\nu_{e})=(0.93\pm 0.30\pm 0.05)\%$
Ablikim:2017omq ; Ablikim:2016rqq . In year 2019, the BESIII collaboration
finished the improved measurements on the branching fractions by using
$e^{+}e^{-}$ annihilation data corresponding to an integrated luminosity of
$3.19~{}{\rm fb}^{-1}$ collected at a center-of-mass energy of $4.178~{}{\rm
GeV}$ and then gave the first determination of $D_{s}\to\eta^{(\prime)}$ TFF,
e.g. ${\cal B}(D_{s}^{+}\to\eta e^{+}\nu_{e})=(2.323\pm 0.063\pm 0.063)\%$,
${\cal B}(D_{s}^{+}\to\eta^{\prime}e^{+}\nu_{e})=(0.824\pm 0.073\pm 0.027)\%$,
$f^{\eta}_{+}(0)|V_{cs}|=0.4455\pm 0.0053\pm 0.0044$, and
$f^{\eta^{\prime}}_{+}(0)|V_{cs}|=0.477\pm 0.049\pm 0.011$ Ablikim:2019rjz .
There are large discrepancies for
$D_{s}^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$ among different experimental
collaborations. With more and more data accumulated in the near future, the
experimental precision shall be greatly improved and the gap among different
measurements could be shrunk.

Theoretically, the decay widths or branching fractions for the semileptonic
decay $D_{s}^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$ depends heavily on the
precision of the $D_{s}\to\eta^{(\prime)}$ TFFs. At present, the
$D_{s}\to\eta^{(\prime)}$ TFFs have been studies under various approaches,
such as the lattice QCD (LQCD) Bali:2014pva , the traditional and covariant
light-front quark model (LFQM) Verma:2011yw ; Cheng:2017pcq ; Wei:2009nc , the
constituent quark model (CQM) Melikhov:2000yu , the covariant confined quark
model (CCQM) Soni:2018adu ; Ivanov:2019nqd , the light-cone sum rules (LCSR)
Duplancic:2015zna ; Offen:2013nma , the QCD sum rules (QCD SR)
Colangelo:2001cv . The LCSR approach is based on the operator product
expansion (OPE) near the light-cone $x^{2}\rightsquigarrow 0$ and
parameterizes all the non-perturbative dynamics into the light-cone
distribution amplitudes (LCDAs), which have been applied for dealing with many
semileptonic decay processes Duplancic:2008ix ; Momeni:2019uag ; Descotes-
Genon:2019bud ; Cheng:2019tgh ; Momeni:2020zrb ; Emmerich:2018rug ;
Momeni:2018tjf ; Shen:2016hyv ; Straub:2015ica ; YanJun:2011rn ; Wang:2010tz ;
Wang:2008bw ; Beneke:2018wjp ; Du:2003ja ; Yang:2005bv ; Singh:2018yvt ;
Fu:2018yin . One may observe that the predicted values of
$D_{s}\to\eta^{(\prime)}$ TFFs behave differently from various groups. Those
discrepancies indicate that it is important to improve the accuracy of
theoretical calculation. In this paper, we will calculate the
$D_{s}\to\eta^{(\prime)}$ TFFs by using the LCSR approach up to next-to-
leading order (NLO) QCD corrections.

It is well known that the states $\eta$ and $\eta^{\prime}$ are considered as
candidates for mixing. The mixing among pseudoscalar mesons is of great
theoretical interests and significance for understanding the dynamics and
hadronic structures, which is caused by the QCD anomaly and related to the
chiral symmetry breaking. Thus, one can gain a better insight into the
dynamics if the mixing parameters are more accurately determined. The
semileptonic decays $D_{s}^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$ probe the
$s\bar{s}$ components of $\eta$ and $\eta^{\prime}$, which can well separate
the strong and weak effects in theory, are expected to be sensitive to
$\eta-\eta^{\prime}$ mixing angle Anisovich:1997dz . Many measurements on the
processes, where $\eta$ and $\eta^{\prime}$ are involved, have been carried
out to fix the mixing parameters. The $\eta-\eta^{\prime}$ mixing can be
described in different forms. To investigate the $\eta-\eta^{\prime}$ mixing,
two schemes have been suggested in the literature Ball:2007hb ; Ball:1995zv ;
Feldmann:1999uf ; Huang:2006as ; Ke:2009mn ; DeFazio:2000my ; Choi:2010zb ,
i.e. the singlet-octet (SO) mixing scheme and the quark-flavor (QF) mixing
scheme. The SO mixing angle $\theta$ between $\eta$ and $\eta^{\prime}$ is
known to be in the range of $\theta\in[-10^{\circ},-23^{\circ}]$. In this
scheme, The $\eta$ and $\eta^{\prime}$ are the mixtures of the flavor SU(3)
octet $\eta_{8}$ and single $\eta_{0}$ states:

$\displaystyle\left(\begin{array}[]{l}\eta\\\
{\eta^{\prime}}\end{array}\right)=\left(\begin{array}[]{l}\cos\theta~{}~{}~{}-\sin\theta\\\
\sin\theta~{}~{}~{}~{}~{}\cos\theta\end{array}\right)\left(\begin{array}[]{l}{\eta_{8}}\\\
{\eta_{0}}\end{array}\right),$ (7)

where $\eta_{8}=(u\bar{u}+d\bar{d}-2s\bar{s})/\sqrt{6}$ and
$\eta_{0}=(u\bar{u}+d\bar{d}+s\bar{s})/\sqrt{3}$. Analogously, information
could be gathered on the mixing scheme in the QF basis, which is consists with
the form $\eta$ and $\eta^{\prime}$ states as combinations of
$|\eta_{q}\rangle=|\bar{u}u+\bar{d}d\rangle/\sqrt{2}$ and
$|\eta_{s}\rangle=|\bar{s}s\rangle$:

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

It has been shown that in this scheme a single angle is essentially required.
In year 2007, the KLOE Collaboration provides the value $\varphi=(41.5\pm
0.3_{\rm stat}\pm 0.7_{\rm syst}\pm 0.6_{\rm th})^{\circ}$ by extracting the
pseudoscalar mixing angle $\varphi$ in the QF basis by measuring the radio
${\cal B}(\phi\to\eta^{\prime}\gamma)/{\cal B}(\phi\to\eta\gamma)$
Ambrosino:2006gk . Some theoretical groups have calculated the single mixing
angle $\varphi$ Ball:2007hb ; Duplancic:2015zna ; Colangelo:2001cv ;
Azizi:2010zj , their predicted values are within the range of
$\varphi\in[39^{\circ},41.8^{\circ}]$. One can put forward the ratio ${\cal
R}_{\eta^{\prime}/\eta}={\cal
B}(D_{s}\to\eta^{\prime}\ell^{+}\nu_{\ell})/{\cal
B}(D_{s}\to\eta\ell^{+}\nu_{\ell})$ to access the $\eta-\eta^{\prime}$ mixing
angle through the ratio of the TFFs
$f_{+}^{\eta}(q^{2})/f_{+}^{\eta^{\prime}}(q^{2})$ Colangelo:2001cv , which
are related to the $\eta-\eta^{\prime}$ mixing scheme. In particular,
information could be gathered on the mixing scheme in the QF basis. In this
paper, we will use the QF basis with the single mixing angle $\varphi$ to
analyze the $D_{s}\to\eta^{(\prime)}$ decay modes, and the corresponding TFFs
satisfy the relation

$\displaystyle\tan\varphi=\frac{|f_{+}^{\eta}(q^{2})|}{|f_{+}^{\eta^{\prime}}(q^{2})|}.$
(9)

One usually takes the large recoil point of the squared momentum transfer
$q^{2}=0$ to do the calculation, i.e.
$\tan\varphi={|f_{+}^{\eta}(0)|}/{|f_{+}^{\eta^{\prime}}(0)|}$. A more
accurate $D_{s}^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$ LCSR analysis is
important.

The $\eta$ and $\eta^{\prime}$ mesons full of rich phenomenology, which are
predominantly flavor-singlet states. This means that their wave functions are
approximately symmetric in the three lightest quark types (up, down, and
strange), which build up the light-hadron spectroscopy. Thus, the LCDAs for
$\eta^{(\prime)}$-meson, as one of the most important parameters, composed by
$s\bar{s}$ are significant to the $D_{s}\to\eta^{(\prime)}$ TFFs, which can be
expanded as a Gegenbauer polynomial series:

$\displaystyle\phi_{2;\eta^{(\prime)}}(u,\mu)=6u\bar{u}\left[1+\sum_{n=1}^{\infty}a^{n}_{2;\eta^{(\prime)}}(\mu)C^{3/2}_{n}(\xi)\right],$
(10)

where $a^{n}_{2;\eta^{(\prime)}}(\mu)$ stands for the $n{\rm th}$-order
Gegenbauer moment, $\bar{u}=(1-u)$ and $\xi=(2u-1)$. When the factorization
scale $\mu$ is large enough, the twist-2 LCDAs
$\phi_{2;\eta^{(\prime)}}(u,\mu)$ tends to the asymptotic form
$\phi_{2;\eta^{(\prime)}}(u,\infty)=6u\bar{u}$. There are some theoretical and
experimental predictions for the Gegenbauer moments $a^{n}_{2;\eta}(\mu)$,
such as the fitting results coming from CLEO collaboration
$a^{2}_{2;\eta}(\mu_{0})=-0.07\pm 0.03$ Gronberg:1997fj , the fitting results
from BABAR result $a^{2}_{2;\eta}(\mu_{0})=-0.05\pm 0.02$ BABAR:2011ad , the
results predicted by Kroll and Passek-Kumericki
$a^{2}_{2;\eta}(\mu_{0})=-0.05\pm 0.02$ Kroll:2013iwa , and the fitting
results from a sum rule analysis $a_{2;\eta}^{2}(\mu_{0})=0.25\pm 0.15$
Offen:2013nma . By taking the approximation with $\pi,K$-meson, Ball and
Zwicky predicted $a^{2}_{2;\eta}(\mu_{0})=0.115$ and
$a^{4}_{2;\eta}(\mu_{0})=-0.015$ Ball:2004ye . At present, few works have been
done to calculate the second and higher order moment. Particularly, the
twist-2 LCDA of $\eta^{\prime}$-meson is rarely studied. So it is important
and meaningful to make a more accurate calculation on the second order moment
and the higher order moments for $\eta^{(\prime)}$-meson LCDAs within the QCD
sum rule approach.

An effective way to calculate the $n$th-order moments of the
$\eta^{(\prime)}$-meson LCDAs is to use the following definition,

$\displaystyle\langle\xi_{2;\eta^{(\prime)}}^{n}\rangle|_{\mu}=\int_{0}^{1}du~{}\xi^{n}~{}\phi_{2;\eta^{(\prime)}}(u,\mu).$
(11)

$\langle\xi_{2;\eta^{(\prime)}}^{n}\rangle|_{\mu}$ can be calculated by using
the Shifman-Vainshtein-Zakharov (SVZ) sum rules Shifman:1978by ; Huang:1986wm
; Govaerts:1983ka ; Huang:1989gv ; Novikov:1983gd ; Hubschmid:1982pa . In
which, the perturbative QCD is established on the assumption that the
perturbative vacuum and the short-distance interaction are not affected by the
long-distance structure of the non-Abelian gauge field. The QCD physical
vacuum contains a series of vacuum condensates, such as the quark condensate
$\langle q\bar{q}\rangle$, the gluon condensate $\langle G^{2}\rangle$, and
etc.. These vacuum condensates reflect the non-perturbation characteristics of
QCD. The QCD sum rules based on the background field theory (BFTSR) method
gives a possible way to consider the non-perturbation effect, which also
provides a systematic description of these vacuum condensates from the field
theory point of view. At present, the BFTSR has been used in calculating the
twist-2 or twist-3 LCDAs for $\pi$, $K$, $D$, $D_{s}$, $a_{0}$, $K_{0}^{*}$,
$f_{0}$, $\rho$, $J/\psi$, $a_{1}(1260)$-mesons Zhong:2014jla ; Fu:2016yzx ;
Fu:2018vap ; Zhong:2014fma ; Zhong:2016kuv ; Zhong:2011jf ; Han:2013zg ;
Huang:2004tp ; Huang:2005av ; Zhong:2011rg ; Zhang:2017rwz ; Zhang:2021wnv ;
Hu:2021lkl . Besides, other methods in studying the LCDAs can be found in
Refs. Cloet:2013tta ; Polyakov:2020cnc ; Cheng:2020vwr ; Wang:2019msf ;
Zhang:2017bzy ; Zuo:2011sk ; Wu:2010zc ; Khodjamirian:2006st ; Zuo:2006re . In
this paper, we will calculate the $\eta^{(\prime)}$-meson twist-2 LCDAs within
the BFTSR approach for the first time.

The rest of the paper are organized as follows. In Sec. II, we present the
basic idea of the QCD background field theory, the detailed BFTSR procedures
for calculating the moments of $\phi_{2;\eta^{(\prime)}}(u,\mu)$, the
branching fractions and the transition form factors involved in the
semileptonic decay $D_{s}^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$. In Sec.
III, we present our numerical results and make a detailed comparison with
other experimental and theoretical predictions. Section IV is reserved for a
summery. The intermediate processes for calculating the moments of
$\eta^{(\prime)}$-meson LCDA within the BFTSR and the basic definition of
$\eta^{(\prime)}$-meson twist-2, 3, 4 LCDAs are given in the Appendixes A, B
and C.

##  II Calculation Technology

###  II.1 Basic idea and formulas for background field theory

The important aspect for the SVZ sum rules approach is the OPE, which has been
introduced by using the QCD physical vacuum $\langle 0|{\cal
O}_{N}|0\rangle_{\rm phys.}$. Its special property is the non-perturbation
effect, which can be described by the classical background field satisfying
the equations of motion. The main idea of the background field theory method
is to describe the non-perturbation effect with the classical background field
satisfying the equation of motion and to describe the quantum fluctuation,
namely the perturbation effect, on this basis within the frame of the quantum
field theory. More specifically, one can use the following substitution in the
theoretical Lagrangian and Green’s functions Makeenko:1979pb ; Shifman:1978zq
; Huang:1986wm

$\displaystyle{\cal A}_{\mu}^{a}(x)\to{\cal
A}_{\mu}^{a}(x)+\phi_{\mu}^{a}(x),$ $\displaystyle\psi(x)\to\psi(x)+\eta(x),$
(12)

where ${\cal A}_{\mu}^{a}(x)$ and $\psi(x)$ stand for the gluon and quark
background fields, $\phi_{\mu}^{a}(x)$ and $\eta(x)$ represent their quantum
fluctuations, respectively. In the presence of background fields, the
quantization of $\phi_{\mu}^{a}(x)$ and $\eta(x)$ has been completed in Ref.
Ambjorn:1982bp . Among them, the gluon quantum field satisfies the background
field gauge,

$\displaystyle D_{\mu}^{ab}({\cal A})\phi_{b}^{\mu}=0,$ $\displaystyle
D_{\mu}^{ab}({\cal A})=\delta^{ab}\partial_{\mu}-g_{s}f^{abc}{\cal
A}_{\mu}^{c},$ (13)

where the color indices $a,b,c=(1,2,...,8)$, $g_{s}$ is the coupling constant
of strong interactions, $f^{abc}$ is the structure constant of the SUf(3)
group. The advantage of choosing the background field gauge is that it makes
the theory be invariant under the background field gauge, making the
calculated physical quantity independent of the gauge. With the help of Eq.
(12), one can obtain the following effective Lagrangian Huang:1989gv

$\displaystyle{\cal L}_{\rm eff}$ $\displaystyle={\cal L}_{\rm QCD}({\cal
A},\psi)+{\cal L}({\rm
ghosts})+\bar{\eta}(i/\\!\\!\\!\\!D-m)\eta+\frac{1}{2}\phi_{\mu}^{a}\bigg{\\{}g^{\mu\nu}D_{ac}^{2}-\Big{(}1-\frac{1}{\alpha}\Big{)}[D^{\mu}{D^{\nu}}]_{ac}$
$\displaystyle+2{g_{s}}f^{abc}G_{b}^{\mu\nu}\bigg{\\}}\phi_{\nu}^{c}~{}+~{}g_{s}(\bar{\psi}/\\!\\!\\!\phi^{a}T^{a}\eta+\bar{\eta}/\\!\\!\\!\phi^{a}T^{a}\eta)~{}-~{}g_{s}^{2}f^{adf}f_{abc}{\cal
A}_{d}^{\mu}\phi_{f}^{\nu}\phi_{\mu}^{b}\phi_{\nu}^{c}-g_{s}f^{abc}$
$\displaystyle\times(\partial_{\mu}\phi_{\nu}^{a})\phi_{b}^{\mu}\phi_{c}^{\nu}-\frac{1}{4}g_{s}^{2}f^{abc}f_{acf}\phi_{b}^{\mu}\phi_{c}^{\nu}\phi_{\mu}^{d}\phi_{\nu}^{f}+{g_{s}}\bar{\eta}{/\\!\\!\\!\phi^{a}}{T^{a}}\eta,~{}$
(14)

where the ${\cal L}_{\rm QCD}({\cal A},\psi)$ with ${\cal A}_{\mu}^{a}(x)$,
$\psi(x)$ has the usual form of QCD Lagrangian and can be minimized to zero
when the classical fields ${\cal A}_{\mu}^{a}(x)$ and $\psi(x)$ are the
solutions of the equation of motion, and $\alpha$ is the gauge-fixing
parameter. The ${\cal L}({\rm ghosts})$ is the contribution of the ghost
particle term. Especially, in order to describe the various quark-antiquark
pairs and gluons in a vacuum, one can follow the classical QCD Lagrangian
Huang:1989gv

$\displaystyle{\cal L}_{\rm QCD}({\cal
A},\psi)=-\frac{1}{4}G_{\mu\nu}^{a}G^{a\mu\nu}+\bar{\psi}(i/\\!\\!\\!\\!D-m)\psi,$
(15)

with $G_{\mu\nu}^{a}=\partial_{\mu}{\cal A}_{\nu}^{a}-\partial_{\nu}{\cal
A}_{\mu}^{a}+g_{s}f^{abc}{\cal A}_{\nu}^{b}{\cal A}_{\mu}^{c}$ stands for the
gluon field strength tensor. The gluon field ${\cal A}_{\mu}^{a}(x)$ and quark
field $\psi(x)$ satisfy the QCD equations of motion,

$\displaystyle(i/\\!\\!\\!\\!D-m)\psi(x)=0,$
$\displaystyle\widetilde{D}_{\mu}^{ab}G^{\nu\mu}_{b}(x)=g_{s}\bar{\psi}(x)\gamma^{\nu}T^{a}\psi(x),$
(16)

where $D_{\mu}=\partial_{\mu}-ig_{s}T^{a}{\cal A}_{\mu}^{a}(x)$ with
$a,b,c=(1,2,...,8)$ and
$\widetilde{D}_{\mu}^{ab}=\delta^{ab}\partial_{\mu}-g_{s}f^{abc}{\cal
A}_{\mu}^{c}(x)$ are fundamental and adjoint representations of the gauge
covariant derivative, respectively. As an advantage of using the background
field theory, one can take different gauges for dealing with the quantum
fluctuations and background fields. More specifically, one can adopt the
background gauge, i.e. $\widetilde{D}^{AB}_{\mu}\phi^{B\mu}(x)=0$ for the
gluon quantum field BG3 ; BG4 , the Schwinger gauge or the fixed-point gauge,
i.e. $x^{\mu}\mathcal{A}^{A}_{\mu}(x)=0$ for the background field
Shifman:1980ui .

As the background field satisfies the motion equation (16), at least one
background field is included in the coupling between the quantum fluctuation
field and the background field. In the effective Lagrangian, there is no
contribution from the vertex
${g_{s}}\bar{\psi}{\gamma^{\mu}}\phi_{\mu}^{a}{T^{a}}\psi$, and it only has
the vertex ${g_{s}}\bar{\psi}/\\!\\!\\!\phi_{\mu}^{a}{T^{a}}\eta$ and its
conjugate contributions. According to the effective Lagrangian, the quantum
quark and gluon propagators in the background field are Huang:2004tp

$\displaystyle S_{F}(x)=i[i\gamma^{\mu}D_{\mu}-m]^{-1},$ (17) $\displaystyle
S_{\mu\nu}^{ab}(x)=i[g_{\mu\nu}({D^{2}})^{ab}+2{g_{s}}f^{abc}G_{\mu\nu}^{c}]^{-1},$
(18)

where the gauge-fixing parameter is taken as $\alpha=1$. Within the framework
of BFT, the quark propagator will be affected by the background quark and/or
gluon fields, which satisfies the equation

$\displaystyle(i/\\!\\!\\!\\!D-m)S_{F}(x,0)=\delta^{4}(x).$ (19)

If taking

$\displaystyle S_{F}(x,0)=(i/\\!\\!\\!\\!D+m)\mathcal{D}(x,0),$ (20)

Eq. (19) can be changed as

$\displaystyle(\Box-\mathcal{P}_{\mu}\partial^{\mu}-\mathcal{Q}+m^{2})\mathcal{D}(x,0)=\delta^{4}(x),$
(21)

where $\Box=\partial^{2}$, and

$\displaystyle\mathcal{P}_{\mu}=2i\mathcal{A}_{\mu}(x),$
$\displaystyle\mathcal{Q}=\gamma^{\nu}\gamma^{\mu}\mathcal{A}_{\nu}(x)\mathcal{A}_{\mu}(x)+i\gamma^{\nu}\gamma^{\mu}\partial_{\nu}\mathcal{A}_{\mu}(x).$
(22)

Moreover, after applying the fixed-point gauge, the gluon background field can
be expressed by using the gauge invariant
$G_{\mu\nu;\alpha_{1}\cdots\alpha_{n}}$ as

$\displaystyle\mathcal{A}_{\mu}(x)$
$\displaystyle=\frac{1}{2}x^{\nu}G_{\nu\mu}+\frac{1}{3}x^{\nu}x^{\alpha}G_{\nu\mu;\alpha}+\frac{1}{8}x^{\nu}x^{\alpha}x^{\beta}G_{\nu\mu;\alpha\beta}+\frac{1}{30}x^{\nu}x^{\alpha}x^{\beta}x^{\gamma}G_{\nu\mu;\alpha\beta\gamma}$
$\displaystyle+\frac{1}{144}x^{\nu}x^{\alpha}x^{\beta}x^{\gamma}x^{\delta}G_{\nu\mu;\alpha\beta\gamma\delta}+\cdots,$
(23)

where $G_{\mu\nu;\alpha_{1}\cdots\alpha_{n}}$ is the notation for
$(D_{\alpha_{1}}\cdots D_{\alpha_{n}}G_{\mu\nu})(0)$, where the indexes
$\alpha_{1}\cdots\alpha_{n}$ indicates the covariant derivative up to $n$-th
order. Substituting Eq. (23) into Eq. (21), we obtain the expressions for
$\mathcal{D}(x,0)$. By further using Eq. (20), we obtain the required quark
propagators in the background field, i.e.,

$S_{F}(x,0)=S_{F}^{0}(x,0)+S_{F}^{2}(x,0)+S_{F}^{3}(x,0)+\sum_{i=1}^{2}S_{F}^{4(i)}(x,0)+\sum_{i=1}^{3}S_{F}^{5(i)}(x,0)+\sum_{i=1}^{5}S_{F}^{6(i)}(x,0).$
(24)

We present the quark propagators with various gauge invariant tensors
$G_{\mu\nu;\alpha_{1}\cdots\alpha_{n}}$ that shall result in up to dimension-
six operators in the Appendix A.1. Because the fixed-point gauge violates the
translation invariance, the quark propagator from $x$ to $0$, $S_{F}(0,x)$,
can not be directly obtained by applying the replacement $x\to-x$ in Eq. (24).
However, it can be related with $S_{F}(x,0)$ via the relation Hubschmid:1982pa

$\displaystyle S_{F}(0,x|\mathcal{A})=CS^{\rm T}_{F}(x,0|-\mathcal{A}^{\rm
T})C^{-1},$ (25)

where $C$ stands for the charge conjugation matrix and the symbol $T$
indicates transposition of both the Dirac and the color matrices.

Furthermore, one will encounter the vertex operator
$\Gamma(z\cdot\tensor{D})^{n}$ with $\Gamma$ indicates all kinds of Dirac
matrices for heavy/light-meson twist-2, 3 LCDAs. Generally, we have the
following expansion Zhong:2014jla

$\displaystyle(z\cdot\tensor{D})^{n}=(z\cdot\overrightarrow{D}-z\cdot\overleftarrow{D})^{n}=(z\cdot\tensor{\partial}+z\cdot
B)^{n}+\cdots,$ (26)

where the ellipsis denotes for the higher-order terms, which are irrelevant
for our present analysis and

$\displaystyle z\cdot B$ $\displaystyle=$ $\displaystyle-2iz\cdot\mathcal{A}$
(27) $\displaystyle=$ $\displaystyle-
ix^{\mu}z^{\nu}G_{\mu\nu}-\frac{2i}{3}x^{\mu}x^{\rho}z^{\nu}G_{\mu\nu;\rho}-\frac{i}{4}x^{\mu}x^{\rho}x^{\sigma}z^{\nu}G_{\mu\nu;\rho\sigma}-\frac{i}{15}x^{\mu}x^{\rho}x^{\sigma}x^{\lambda}z^{\nu}G_{\mu\nu;\rho\sigma\lambda}$
$\displaystyle-$
$\displaystyle\frac{i}{72}x^{\mu}x^{\rho}x^{\sigma}x^{\lambda}x^{\tau}z^{\nu}G_{\mu\nu;\rho\sigma\lambda\tau}+\cdots.$

We can expand the operator $(z\cdot\tensor{D})^{n}$ into series of the
operators $(z\cdot\tensor{\partial})^{n}$ and $G_{\mu\nu;\rho\cdots}$. For the
purpose, we first expand $(z\cdot\tensor{D})^{n}$ as

$\displaystyle(z\cdot\tensor{D})^{0}$ $\displaystyle=$ $\displaystyle 1,$
$\displaystyle(z\cdot\tensor{D})^{1}$ $\displaystyle=$ $\displaystyle
z\cdot\tensor{\partial}+z\cdot B,$ $\displaystyle(z\cdot\tensor{D})^{2}$
$\displaystyle=$
$\displaystyle(z\cdot\tensor{\partial})^{2}+2(z\cdot\tensor{\partial})(z\cdot\underline{B})+(z\cdot
B)^{2},$ $\displaystyle(z\cdot\tensor{D})^{3}$ $\displaystyle=$
$\displaystyle(z\cdot\tensor{\partial})^{3}+3(z\cdot\tensor{\partial})^{2}(z\cdot\underline{B})+\left[(z\cdot\partial)^{2}(z\cdot
B)\right]+3(z\cdot\tensor{\partial})(z\cdot\underline{B})^{2}+(z\cdot B)^{3},$
(28) $\displaystyle\cdots\cdots,$

where, the operator “underline” below symbol “$B$” (or “$x$” in the Appendix
A.2) indicates that the operation $\tensor{\partial}$ does not act on it. In
deriving those equations, the following equation has been adopted,

$\displaystyle(z\cdot\tensor{\partial})^{n}(z\cdot
B)=\sum^{n}_{k=0}\frac{n!}{k!(n-k!)}(z\cdot\tensor{\partial})^{n-k}\left[(z\cdot\partial)^{k}(z\cdot\underline{B})\right].$
(29)

By keeping only those terms that shall leads to operators up to dimension-six,
we can obtain the full expression for the vertex operator, which are listed in
the Appendix A.2.

###  II.2 SVZ sum rules for the moments of $\phi_{2;\eta^{(\prime)}}(x,\mu)$

Following the traditional way of constructing the light pseudoscalar meson, we
take the following definition

$\displaystyle\langle 0|\bar{s}(0){\cal
C}_{s}/\\!\\!\\!z\gamma_{5}[z,-z]{(iz\cdot\tensor{D})^{n}}s(0)|\eta^{(\prime)}(q)\rangle=i{(z\cdot
q)}^{n+1}f_{\eta^{(\prime)}}\langle\xi^{n}_{2;\eta^{(\prime)}}\rangle|_{\mu},$
(30)

where $z_{\mu}$ stand for the light-like vector and $[z,-z]$ is the path-
ordered gauge connection, $f_{\eta^{(\prime)}}$ are the $\eta$ and
$\eta^{\prime}$-meson decay constant,
$(iz\cdot\tensor{D})^{n}=(iz\cdot\overrightarrow{D}-iz\cdot\overleftarrow{D})^{n}$.
Basis on the QF scheme, the flavour content is ${\cal C}_{s}=({\cal
C}_{1}-\sqrt{2}{\cal C}_{8})/\sqrt{3}$ with SO basis ${\cal
C}_{1}=\mathbf{1}/\sqrt{n_{f}}$ and ${\cal C}_{8}=\lambda_{8}/\sqrt{2}$
Duplancic:2015zna ; Ball:2007hb . In which the $\lambda_{i}$ is the standard
SUf(3) Gell-Mann matrix and $\mathbf{1}$ is the $3\times 3$ unit matrix. As a
special case, the 0th-order LCDA’s moment for Eq. (11) satisfies the
normalization condition

$\displaystyle\langle\xi^{0}_{2;\eta^{(\prime)}}\rangle|_{\mu}=\int_{0}^{1}du\phi_{2;\eta^{(\prime)}}(u,\mu)=1.$
(31)

To derive the SVZ sum rules for the moments
$\langle\xi^{n}_{2;\eta^{(\prime)}}\rangle|_{\mu}$, we introduce the following
correlation function,

$\displaystyle\Pi_{2;\eta^{(\prime)}}^{(n,0)}(z,q)$ $\displaystyle=i\int
d^{4}xe^{iq\cdot x}\langle 0|T\\{J_{n}(x),J_{0}^{\dagger}(0)\\}|0\rangle$
$\displaystyle=(z\cdot q)^{n+2}I_{2;\eta^{(\prime)}}^{(n,0)}(q^{2}),$ (32)

where $J_{n}(x)=\bar{s}(x){\cal
C}_{s}/\\!\\!\\!z\gamma_{5}(iz\cdot\tensor{D})^{n}s(x)$,
$J_{0}^{\dagger}(0)=\bar{s}(0){\cal C}_{s}/\\!\\!\\!z\gamma_{5}s(0)$ and
$z^{2}\rightsquigarrow 0$. Only even moments are non-zero and the odd moments
of the LCDA are zero because of G-parity, i.e. $n=(0,2,4,6,...)$ will
contribute to the final results 111This point can also been seen in Eq. (16)
of Ref. Ball:2007hb .. At one hand, in deep Euclidean region $q^{2}\ll 0$, one
can apply the OPE for the correlation function Eq. (32).

$\displaystyle\Pi_{2;\eta^{(\prime)}}^{(n,0)}(z,q)=i\int d^{4}xe^{iq\cdot
x}{\rm Tr}[{\cal C}_{s}{\cal C}_{s}]\Big{\\{}\\!$ $\displaystyle-{\rm
Tr}\langle
0|S_{F}^{s}(0,x)/\\!\\!\\!z{{\gamma}_{5}}(iz\cdot\tensor{D})^{n}S_{F}^{s}(x,0)/\\!\\!\\!z\gamma_{5}|0\rangle$
$\displaystyle+{\rm Tr}\langle
0|\bar{s}(x)s(0)/\\!\\!\\!z{{\gamma}_{5}}{(iz\cdot\tensor{D})^{n}}S_{F}^{s}(x,0)/\\!\\!\\!z{{\gamma}_{5}}|0\rangle$
$\displaystyle+{\rm Tr}\langle
0|S_{F}^{s}(0,x)/\\!\\!\\!z{{\gamma}_{5}}{(iz\cdot\tensor{D})^{n}}\bar{s}(0)s(x)/\\!\\!\\!z{{\gamma}_{5}}|0\rangle$
$\displaystyle+\cdots\Big{\\}},$ (33)

where ${\rm Tr}[{\cal C}_{s}{\cal C}_{s}]=1$. In the detailed OPE calculation,
we adopt the $\overline{\rm MS}$-scheme to deal with the infrared divergences.
Lorentz invariant scalar function $I_{2;\eta^{(\prime)}}^{(n,0)}(q^{2})$ in
Eq. (32) depends on the condensation parameter and will encounter the vacuum
matrix elements of the following form,

$\displaystyle\langle 0|G_{\mu\nu}^{a}G_{\mu\nu}^{b}|0\rangle,$
$\displaystyle\langle
0|G_{\mu\nu}^{a}G_{\rho\sigma}^{b}G_{\lambda\tau}^{c}|0\rangle,$
$\displaystyle\langle
0|G_{\mu\nu;\lambda}^{a}G_{\rho\sigma;\tau}^{b}|0\rangle,$
$\displaystyle\langle
0|G_{\mu\nu}^{a}G_{\rho\sigma;\lambda\tau}^{b}|0\rangle,$
$\displaystyle\langle
0|G_{\mu\nu;\lambda\tau}^{a}G_{\rho\sigma}^{c}|0\rangle,$
$\displaystyle\langle 0|\bar{q}_{\alpha}^{a}(x)q_{\beta}^{b}(y)|0\rangle,$
$\displaystyle\langle
0|\bar{q}_{\alpha}^{a}(x)q_{\beta}^{b}(y)G_{\mu\nu}^{A}|0\rangle,$
$\displaystyle\langle
0|\bar{q}_{\alpha}^{a}(0)q_{\beta}^{b}(0)G_{\mu\nu;\rho}^{A}|0\rangle.$ (34)

The full expression for the vacuum condensates which one may encounter in the
following calculations are listed in the Appendix A.3, which can also be found
in our previous work Zhong:2014jla . On the other side, the correlation
function (32) can be treated by inserting a complete set of intermediate
hadronic states in physical region to obtain its hadronic representation

$\displaystyle{\rm Im}I_{2;\eta^{(\prime)},{\rm Had}}^{(n,0)}(q^{2})$
$\displaystyle=\pi\delta(q^{2}-m_{\eta^{(\prime)}}^{2})f_{\eta^{(\prime)}}^{2}\langle\xi^{n}_{2;\eta^{(\prime)}}\rangle|_{\mu}+\pi\frac{3}{4\pi^{2}(n+1)(n+3)}\theta(q^{2}-s_{\eta^{(\prime)}}),$
(35)

in which $s_{\eta^{(\prime)}}$ stand for the continuum threshold for the
lowest continuum state. The first term is the contribution of
$\eta^{(\prime)}$-meson poles, and the second term is the contribution of
continuum states above poles. Then both of the OPE part and the hadronic
representation of the invariant function
$I_{2;\eta^{(\prime)}}^{(n,0)}(q^{2})$ can be marched with the following
dispersion relation

$\displaystyle\frac{1}{\pi}\int_{4m_{s}^{2}}^{\infty}{ds}\frac{{\rm
Im}I_{2;\eta^{(\prime)},{\rm
Had}}^{(n,0)}(s)}{s-q^{2}}=I_{2;\eta^{(\prime)},{\rm
QCD}}^{(n,0)}({{q}^{2}}).$ (36)

The Borel transform helps to reduce the contribution from the continuum states
on the left of Eq. (36) and the contribution of the higher-dimension
condensates on the right, and finally the sum rules can be obtained,

$\displaystyle\frac{1}{\pi}\frac{1}{M^{2}}\int dse^{-s/M^{2}}{\rm
Im}I_{2;\eta^{(\prime)},{\rm had}}^{(n,0)}(s)=\hat{\cal
B}_{M^{2}}I_{2;\eta^{(\prime)},{\rm QCD}}^{(n,0)}(q^{2}),$ (37)

where $M^{2}$ is Borel parameter, $\hat{\cal B}_{M^{2}}$ is Borel
transformation operator,

$\displaystyle\hat{\cal
B}_{M^{2}}=\mathop{\lim}\limits_{\scriptstyle-q^{2},n\to\infty\hfill\atop\scriptstyle{-q^{2}/n}={M^{2}}\hfill}\frac{1}{(n-1)!}(-q^{2})^{n}\Bigg{(}-\frac{d}{d(-q^{2})}\Bigg{)}^{n}.$
(38)

In order to deal with the $s$-quark mass ($m_{s}$) contribution to
$\langle\xi_{2;\eta^{(\prime)}}^{n}\rangle|_{\mu}$, we take the expansion
according to different orders of $m_{s}^{k}$ with $k=(0,2,4,...)$, i.e.
$I_{m_{s}^{k}}(n,M^{2})$, since $m_{s}$ is closer to $\Lambda_{\rm QCD}$,
which is different from our previous treatment for the heavy quark such as
$q=(c,b)$ in Refs. Fu:2018vap ; Zhong:2014fma . Here, we take the first two
orders of the squared $s$-quark mass, i.e. $I_{m_{s}^{0}}(n,M^{2})$ and
$I_{m_{s}^{2}}(n,M^{2})$. The reason lies in the suppression of
$m_{s}^{4}<0.1\%$ for the third-order, which are quite small and can be safely
neglected. Recently, we have suggested a new method for renormalization of
various moments, i.e. the 0th-order moment
$\langle\xi_{2;\eta^{(\prime)}}^{0}\rangle|_{\mu}$ should be considered to the
total results. The reason lies in the accuracy is often up to dimension-six
condensates and the NLO QCD corrections for the perturbative part instead of
the infinite dimension or infinite-order perturbative parts Zhong:2021epq .
The final expression for $\langle\xi_{2;\eta^{(\prime)}}^{n}\rangle|_{\mu}$
can then be written as

$\displaystyle\frac{f_{\eta^{(\prime)}}^{2}\langle\xi^{n}_{2;\eta^{(\prime)}}\rangle|_{\mu}\langle\xi^{0}_{2;\eta^{(\prime)}}\rangle|_{\mu}}{M^{2}e^{m_{\eta^{(\prime)}}^{2}/M^{2}}}=\frac{1}{\pi}\frac{1}{M^{2}}\int_{4m_{s}^{2}}^{s_{\eta^{(\prime)}}}dse^{-s/M^{2}}\frac{3v^{n+1}}{8\pi(n+1)(n+3)}\left(1+\frac{\alpha_{s}}{\pi}A^{\prime}_{n}\right)\bigg{\\{}[1+(-1)^{n}]$
$\displaystyle\hskip
51.21504pt\times(n+1)\frac{1-v^{2}}{2}+[1+(-1)^{n}]\bigg{\\}}+\frac{2m_{s}\langle\bar{s}s\rangle}{M^{4}}+\frac{\langle\alpha_{s}G^{2}\rangle}{12\pi
M^{4}}~{}\frac{1+n\theta(n-2)}{n+1}\,-\,\frac{8n+1}{9}$ $\displaystyle\hskip
51.21504pt\times\frac{m_{s}\langle g_{s}\bar{s}\sigma
TGs\rangle}{M^{6}}+\frac{\langle
g_{s}\bar{s}s\rangle}{81M^{6}}4(2n+1)-\frac{\langle
g_{s}^{3}fG^{3}\rangle}{48\pi^{2}M^{6}}n\theta(n-2)+\frac{\sum\langle
g_{s}^{2}\bar{q}q\rangle^{2}}{486\pi^{2}M^{6}}\bigg{\\{}\\!\\!-2(51n$
$\displaystyle\hskip
51.21504pt+25)\bigg{(}-\ln\frac{M^{2}}{\mu^{2}}\bigg{)}+3\,(17n+35)+\theta(n-2)\bigg{[}\,2n\bigg{(}-\ln\frac{M^{2}}{\mu^{2}}\bigg{)}-25\,(2n+1)\,\tilde{\psi}(n)$
$\displaystyle\hskip
51.21504pt+\frac{1}{n}(49n^{2}+100n+56)\bigg{]}\bigg{\\}}+I_{m_{s}^{2}}(n,M^{2}).$
(39)

Due to the mass of $s$ quark is heavier than $u$, $d$-quark, the
$I_{m_{s}^{2}}(n,M^{2})$-terms should be considered in this paper, which are

$\displaystyle
I_{m_{s}^{2}}(n,M^{2})=m_{s}^{2}\bigg{\\{}-\frac{\langle\alpha_{s}G^{2}\rangle}{6\pi
M^{6}}\left[\theta(n-2)(n\tilde{\psi}(n)-2)+2n\left(-\ln\frac{M^{2}}{\mu^{2}}\right)-n-2\right]+\frac{\langle
g_{s}^{3}fG^{3}\rangle}{288\pi^{2}M^{8}}$ $\displaystyle\hskip
17.07182pt\times\bigg{\\{}\\!-10\delta^{n0}+\theta(n-2)\left[4n(2n-1)\left(-\ln\frac{M^{2}}{\mu^{2}}\right)-4n\tilde{\psi}(n)+8({n^{2}}-n+1)\right]+\theta(n-4)$
$\displaystyle\hskip
17.07182pt\times[2n(8n-1)\tilde{\psi}(n)-(19{n^{2}}+19n+6)]+8n\,(3n-1)\left(-\ln\frac{M^{2}}{\mu^{2}}\right)-(21n^{2}+53n-6)\bigg{\\}}$
$\displaystyle\hskip 17.07182pt-\frac{\sum\langle
g_{s}^{2}q\bar{q}\rangle^{2}}{972\pi^{2}M^{8}}~{}\bigg{\\{}6{\delta^{n0}}\left[{16\left(-\ln\frac{M^{2}}{\mu^{2}}\right)-3}\right]+\theta(n-2)\bigg{[}8(n^{2}+12n-12)\left(-\ln\frac{M^{2}}{\mu^{2}}\right)-2$
$\displaystyle\hskip
17.07182pt\times(29n+22)\,\tilde{\psi}(n)~{}+~{}4\left(5{n^{2}}-2n-33+\frac{46}{n}\right)\bigg{]}+\theta(n-4)~{}\bigg{[}2\left(56{n^{2}}-25n+24\right)\tilde{\psi}(n)$
$\displaystyle\hskip
17.07182pt-\,\left(139n^{2}+91n+54\right)\bigg{]}~{}+~{}8\left(27n^{2}-15n-11\right)~{}\left(-\ln\frac{M^{2}}{\mu^{2}}\right)-3(63n^{2}+159n-50)\bigg{\\}}$
$\displaystyle\hskip
17.07182pt+\frac{4(n-1)}{3}\frac{m_{s}\langle\bar{s}s\rangle}{M^{6}}+\frac{8n-3}{9}\frac{m_{s}\langle
g_{s}\bar{s}\sigma TGs\rangle}{M^{8}}-\frac{4(2n+1)}{81}\frac{\langle
g_{s}\bar{s}s\rangle^{2}}{M^{8}}\bigg{\\}}.$ (40)

For convenience, we put the detailed terms contribute to the BFTSR
$\langle\xi_{2;\eta^{(\prime)}}^{n}\rangle|_{\mu}$ in the Appendix B. By
taking the index $n$ to zero, we get the sum rule of the 0th-order moment,
which takes the following form,

$\displaystyle\frac{(\langle\xi^{0}_{2;\eta^{(\prime)}}\rangle|_{\mu})^{2}f_{\eta^{(\prime)}}^{2}}{M^{2}e^{m_{\eta^{(\prime)}}^{2}/M^{2}}}=\frac{1}{\pi}\frac{1}{M^{2}}\int_{4m_{s}^{2}}^{s_{\eta^{(\prime)}}}dse^{-s/M^{2}}\frac{v(3-v^{2})}{8\pi}+\frac{2m_{s}\langle\bar{s}s\rangle}{M^{4}}-\frac{m_{s}\langle
g_{s}\bar{s}\sigma
TGs\rangle}{9M^{6}}+\frac{\langle\alpha_{s}G^{2}\rangle}{12\pi M^{4}}$
$\displaystyle\hskip 29.87547pt+\frac{4\langle
g_{s}s\bar{s}\rangle^{2}}{81M^{6}}+\frac{\sum\langle
g_{s}^{2}\bar{q}q\rangle^{2}}{M^{6}}\frac{1}{486{\pi^{2}}}\bigg{[}-50\bigg{(}-\ln\frac{M^{2}}{\mu^{2}}\bigg{)}+105\bigg{]}+m_{s}^{2}\bigg{\\{}\frac{\langle\alpha_{s}G^{2}\rangle}{3\pi
M^{6}}-\frac{\langle g_{s}^{3}fG^{3}\rangle}{72\pi^{2}M^{8}}$
$\displaystyle\hskip 29.87547pt-\frac{\sum\langle
g_{s}^{2}\bar{q}q\rangle^{2}}{972\pi^{2}M^{8}}~{}\Bigg{[}8\left(-\ln\frac{M^{2}}{\mu^{2}}\right)-132\Bigg{]}-\frac{4}{3}\frac{m_{s}\langle\bar{s}s\rangle}{M^{6}}\,-\,\frac{1}{3}\frac{m_{s}\langle
g_{s}\bar{s}\sigma TGs\rangle}{M^{8}}\,-\,\frac{4}{81}\frac{\langle
g_{s}\bar{s}s\rangle^{2}}{M^{8}}\bigg{\\}},$ (41)

where

$\displaystyle\widetilde{\psi}(n)=\psi\bigg{(}\frac{n+1}{2}\bigg{)}-\psi\bigg{(}\frac{n}{2}\bigg{)}+\ln
4,$ (42)

with $v^{2}=1-4m_{s}^{2}/s$ and $A^{\prime}_{0}=0$, $A^{\prime}_{2}=5/3$,
$A^{\prime}_{4}=59/27$, $A^{\prime}_{6}=353/135$ are the NLO correction to the
perturbative part Ball:1996tb . The 0th-order derivative of the digamma
function $\psi(n+1)=\Sigma_{k=1}^{n}1/k-\gamma_{E}$, where the Euler’s
constant $\gamma_{E}=0.557216$. Furthermore, in order to get the Gegenbauer
moment $a_{2;\eta^{(\prime)}}^{n}$, one can expand
$\phi_{2;\eta^{(\prime)}}^{n}(u,\mu)$ into a Gegenbauer polynomial series by
using Eq. (10) and the basic definition of
$\langle\xi^{n}_{2;\eta^{(\prime)}}\rangle|_{\mu}$ (11). Then, one can get the
following relations up to $6$th-order,

$\displaystyle\langle\xi^{2}_{2;\eta^{(\prime)}}\rangle|_{\mu}=\frac{1}{5}+\frac{12}{35}a_{2;\eta^{(\prime)}}^{2}(\mu),$
$\displaystyle\langle\xi^{4}_{2;\eta^{(\prime)}}\rangle|_{\mu}=\frac{3}{35}+\frac{8}{35}a_{2;\eta^{(\prime)}}^{2}(\mu)+\frac{8}{77}a_{2;\eta^{(\prime)}}^{4}(\mu),$
$\displaystyle\langle\xi^{6}_{2;\eta^{(\prime)}}\rangle|_{\mu}=\frac{1}{21}+\frac{12}{77}a_{2;\eta^{(\prime)}}^{2}(\mu)+\frac{120}{1001}a_{2;\eta^{(\prime)}}^{4}(\mu)+\frac{64}{2145}a_{2;\eta^{(\prime)}}^{6}(\mu).$
(43)

Following this method, one can get higher-order Gegenbauer moments.

###  II.3 The semileptonic decay
$D_{s}^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$

All the following calculations are performed under the Standard Model (SM). In
order to derive the full analytical LCSR expressions for the TFFs, we use the
traditional current method to calculate the TFFs. The correlation function for
the TFFs of $D_{s}\to\eta^{(\prime)}$ is defined as Duplancic:2008ix :

$\displaystyle\Pi_{\mu}(p,q)$ $\displaystyle=i\int
d^{4}xe^{iqx}\langle\eta^{(\prime)}(p)|T\\{\bar{s}(x)\gamma_{\mu}c(x),\bar{c}(0)i\gamma_{5}s(0)\\}|0\rangle$
$\displaystyle=\Pi[q^{2},(p+q)^{2}]p_{\mu}+\tilde{\Pi}[q^{2},(p+q)^{2}]q_{\mu}.$
(44)

Following the basic procedure of LCSR approach, the correlation function can
be treated by inserting complete intermediate states with the same quantum
numbers as the current operator $(\bar{c}i\gamma_{5}s)$ in the time-like
${(p+q)}^{2}$-region. After isolating the pole term of the lowest pseudoscalar
$D_{s}$-meson, one can reach the following expression,

$\displaystyle\Pi_{\mu}^{\rm had}(p,q)$
$\displaystyle=\frac{\langle\eta^{(\prime)}(p)|\bar{s}\gamma_{\mu}c|D_{s}(p+q)\rangle\langle
D_{s}(p+q)|\bar{c}i\gamma_{5}q|0\rangle}{m_{D_{s}}^{2}-(p+q)^{2}}$
$\displaystyle+\sum\limits_{\rm
H}{\frac{\langle\eta^{(\prime)}(p)|\bar{s}\gamma_{\mu}c|D_{s}^{\rm
H}(p+q)\rangle\langle D_{s}^{\rm
H}(p+q)|\bar{c}i\gamma_{5}q|0\rangle}{m_{D_{s}^{\rm H}}^{2}-(p+q)^{2}}}$
$\displaystyle=\Pi^{\rm had}[q^{2},(p+q)^{2}]p_{\mu}+\widetilde{\Pi}^{\rm
had}[q^{2},(p+q)^{2}]q_{\mu},$ (45)

with the superscript “had” and “H” stand for the hadronic expression of the
correlation function and higher-excited state of $D_{s}$-meson, respectively.
The decay constant of $D_{s}$-meson can be defined via the relation, $\langle
D_{s}|\bar{c}i\gamma_{5}q|0\rangle=m_{D_{s}}^{2}f_{D_{s}}/m_{c}$. The
definition of transition matrix element for $D_{s}\to\eta^{(\prime)}$ has the
following form

$\displaystyle\langle\eta^{(\prime)}(p)|\bar{s}\gamma_{\mu}c|D_{s}(p+q)\rangle=2f^{\eta^{(\prime)}}_{+}(q^{2})p_{\mu}+\tilde{f}^{\eta^{(\prime)}}(q^{2})q_{\mu}.$
(46)

with
$\tilde{f}^{\eta^{(\prime)}}=f^{\eta^{(\prime)}}_{+}(q^{2})+f^{\eta^{(\prime)}}_{-}(q^{2})$.
Then one can take the imaginary part of the invariant amplitude $\Pi^{\rm
had}[q^{2},(p+q)^{2}]$ and $\widetilde{\Pi}^{\rm had}[q^{2},(p+q)^{2}]$ which
has the following form,

$\displaystyle{\rm Im}\Pi^{\rm
had}(q^{2},s)=\pi\delta(s-m_{D_{s}^{2}})\frac{2m_{D_{s}^{2}}f_{D_{s}}f_{+}^{\eta^{(\prime)}}}{m_{c}}+\pi\rho^{\rm
H}(s)\theta(s-s_{0})$ $\displaystyle{\rm Im}\widetilde{\Pi}^{\rm
had}(q^{2},s)=\pi\delta(s-m_{D_{s}^{2}})\frac{m_{D_{s}^{2}}f_{D_{s}}\tilde{f}_{+}^{\eta^{(\prime)}}}{m_{c}}+\pi\tilde{\rho}^{\rm
H}(s)\theta(s-\tilde{s}_{0})$ (47)

where $\rho^{\rm H}(s)$ and $\tilde{\rho}^{\rm H}(s)$ denote the spectral
density of higher resonance and the continuum states $D_{s}^{\rm H}$, which
can be approximated by invoking the so-called quark-hadron duality ansatz
$\varrho^{\rm H}(s)=\varrho^{\rm QCD}(s)$ with
$\varrho(s)=(\rho(s),\tilde{\rho}(s))$ and the perturbative spectral density
$\rho^{\rm QCD}$ and $\tilde{\rho}^{\rm QCD}$, the usual step-function
$\theta(x)$. Here, the $s_{0}$ ($\tilde{s}_{0}$) are effective parameters,
which characterize the lower limit of continuum state, namely the continuum
threshold parameters. After extracting the pole terms, the
$s_{0}(\tilde{s}_{0})$ can separate the ground state $D_{s}$ and the excited
state $D_{s}^{\rm H}$. When taking the limit for the interval between the two
adjacent excited states, the sum of excited states is transformed into the
integral of continuum states and $s_{0}(\tilde{s}_{0})$ will be changed into
the lower limit of the integration. Traditionally, the continuum threshold are
often taken as the magnitude that is close to the squared mass of the first
excited state. The ground state of $D_{s}$-meson is calculated by the LCSR
approach. The excited states’ contribution can be highly suppressed when
making the Borel transformation. Meanwhile, the continuum states, contribution
is usually required to be less than 30% so as to make the LCSR calculation
more accurate and reliable. Then, one can use a general dispersion relation in
the momentum squared $(p+q)^{2}$, which can establish a relationship with the
QCD parts, i.e.

$\displaystyle\Pi^{\rm QCD}(q^{2},s)$
$\displaystyle=\frac{1}{\pi}\int_{m_{c}^{2}}^{\infty}\frac{{\rm Im}\Pi^{\rm
had}(q^{2},s)}{s-(p+q)^{2}}ds,$ $\displaystyle\widetilde{\Pi}^{\rm
QCD}(q^{2},s)$ $\displaystyle=\frac{1}{\pi}\int_{m_{c}^{2}}^{\infty}\frac{{\rm
Im}\tilde{\Pi}^{\rm had}(q^{2},s)}{s-(p+q)^{2}}ds.$ (48)

Here, we shall only deal with ${\rm Im}\Pi^{\rm had}[q^{2},s]$ for the TFFs
$f^{\eta^{(\prime)}}_{+}(q^{2})$, which are the only TFFs contribute to the
required branching fractions.

On the other hand, when the correlation function (44) is dominated by the
light-like distances, it can be expanded around the light-cone. The light-cone
expansion is performed by integrating out the transverse and minus degrees of
freedom and leaving only the longitudinal momenta of the partons as the
relevant degrees of freedom. The integration over the transverse momenta is
done up to a cutoff, $\mu_{\rm IR}$, all momenta below which are included in
the $\eta^{(\prime)}$-meson LCDAs. Large transverse momenta are calculated in
perturbative theory. Thus, the correlation function can be separated into
perturbative and nonperturbative contributions, both of which depend on the
longitudinal parton momenta and the factorization scale $\mu_{\rm IR}$
Ball:2007hb .

In order to make our result more accurate, we consider both the leading-order
(LO) for all the LCDAs’ part and gluon radiative corrections to the dominant
twist-2, 3 parts of the correlation function. The OPE result for the
correlation function $\Pi^{\rm OPE}$ is then represented as a sum of LO and
NLO parts,

$\displaystyle\Pi^{\rm
OPE}[q^{2},(p+q)^{2}]=F_{0}(q^{2},(p+q)^{2})+\frac{\alpha_{s}C_{F}}{4\pi}F_{1}(q^{2},(p+q)^{2}).$
(49)

To calculate the invariant amplitude $F_{0}(q^{2},(p+q)^{2})$ and
$F_{1}(q^{2},(p+q)^{2})$, one needs to know the expression for the $c$-quark
propagator, i.e.

$\displaystyle\langle 0|T\\{c(x)\bar{c}(0)\\}|0\rangle$
$\displaystyle=i\int\frac{d^{4}k}{(2\pi)^{4}}{e^{-ik\cdot
x}}\frac{m_{c}+/\\!\\!\\!k}{k^{2}-m_{c}^{2}}-ig\int\frac{d^{4}k}{(2\pi)^{4}}e^{-ik\cdot
x}\int_{0}^{1}dv\Bigg{[}\frac{m_{c}+/\\!\\!\\!k}{2(m_{c}^{2}-k^{2})^{2}}$
$\displaystyle\times
G^{\mu\nu}(vx)\sigma_{\mu\nu}+\frac{v}{m_{c}^{2}-k^{2}}x_{\mu}G^{\mu\nu}(vx)\gamma_{\nu}\Bigg{]}.$
(50)

To do the calculation, the expression of twist-2, 3, 4 LCDAs matrix elements
are needed, which are displayed in the Appendix C. To get the final LCSR
expression, we need to use the Borel transformation to transform the variable
$(p+q)^{2}$ into Borel parameter $M^{2}$. Then the expression of the
$D_{s}\to\eta^{(\prime)}$ TFFs up to NLO gluon radiation correction to the
twist-2, 3 LCDAs can be obtained by equating the two types of representation
of the correlation function and by subtracting the contribution from higher
resonances and continuum states, i.e.,

$\displaystyle f^{\eta^{(\prime)}}_{+}(q^{2})$
$\displaystyle=\frac{e^{m_{D_{s}}^{2}/M^{2}}}{2m_{D_{s}}^{2}f_{D_{s}}}\bigg{[}F_{0}(q^{2},M^{2},s_{0})+\frac{\alpha_{s}C_{F}}{4\pi}F_{1}(q^{2},M^{2},s_{0})\bigg{]}.$
(51)

where $F_{0(1)}(q^{2},M^{2},s_{0})$ originates from the OPE result for the LO
(NLO) invariant amplitude $F_{0(1)}(q^{2},(p+q)^{2})$. Finally, the LCSR for
$D_{s}\to\eta^{(\prime)}$ TFFs have the following form,

$\displaystyle
f^{\eta^{(\prime)}}_{+}(q^{2})=\frac{m_{c}^{2}f_{\eta^{(\prime)}}}{2m_{D_{s}}^{2}f_{D_{s}}}e^{m_{D_{s}}^{2}/M^{2}}\\!\\!\int_{u_{0}}^{1}\\!due^{-s(u)/M^{2}}\bigg{\\{}\frac{\phi_{2;\eta^{(\prime)}}(u)}{u}\\!+\frac{1}{2m_{s}m_{c}}\bigg{[}\phi_{3;\eta^{(\prime)}}^{p}(u)\\!+\\!\frac{1}{6}\bigg{(}2\frac{\phi_{3;\eta^{(\prime)}}^{\sigma}(u)}{u}$
$\displaystyle\hskip
7.96674pt-\frac{m_{c}^{2}+q^{2}-u^{2}m_{\eta^{(\prime)}}^{2}}{m_{c}^{2}-q^{2}+u^{2}m_{\eta^{(\prime)}}^{2}}\frac{d}{du}\phi_{3;\eta^{(\prime)}}^{\sigma}(u)+\frac{4um_{\eta^{(\prime)}}^{2}m_{c}^{2}}{(m_{c}^{2}-q^{2}+u^{2}m_{\eta^{(\prime)}}^{2})^{2}}\phi_{3;\eta^{(\prime)}}^{\sigma}(u)\bigg{)}\bigg{]}\\!+\\!\frac{1}{m_{c}^{2}-q^{2}+u^{2}m_{\eta^{(\prime)}}^{2}}$
$\displaystyle\hskip
7.96674pt\times\bigg{[}u\psi_{4;\eta^{(\prime)}}(u)+\left(1-\frac{{2{u^{2}}m_{\eta^{(\prime)}}^{2}}}{m_{c}^{2}-q^{2}+u^{2}m_{\eta^{(\prime)}}^{2}}\right)\int_{0}^{u}dv\psi_{4;\eta^{(\prime)}}(v)-\frac{m_{c}^{2}}{4}\frac{u}{m_{c}^{2}-q^{2}+u^{2}m_{\eta^{(\prime)}}^{2}}\,\bigg{(}\frac{d^{2}}{du^{2}}$
$\displaystyle\hskip
7.96674pt-\frac{6um_{\eta^{(\prime)}}^{2}}{m_{c}^{2}-q^{2}+u^{2}m_{\eta^{(\prime)}}^{2}}\frac{d}{du}\,+\frac{12um_{\eta^{(\prime)}}^{4}}{(m_{c}^{2}-q^{2}+u^{2}m_{\eta^{(\prime)}}^{2})^{2}}\bigg{)}\phi_{4;\eta^{(\prime)}}(u)-\bigg{(}\frac{d}{du}-\frac{2um_{\eta^{(\prime)}}^{2}}{m_{c}^{2}-q^{2}+u^{2}m_{\eta^{(\prime)}}^{2}}\bigg{)}$
$\displaystyle\hskip
7.96674pt\times\bigg{(}\frac{f_{3\eta^{(\prime)}}}{f_{\eta^{(\prime)}}m_{c}}~{}I_{3;\eta^{(\prime)}}(u)~{}+I_{4;\eta^{(\prime)}}(u)\bigg{)}-\frac{{2{u}m_{\eta^{(\prime)}}^{2}}}{m_{c}^{2}-q^{2}+u^{2}m_{\eta^{(\prime)}}^{2}}~{}\bigg{(}u\frac{d}{du}-\frac{{2{u^{2}}m_{\eta^{(\prime)}}^{2}}}{m_{c}^{2}-q^{2}+u^{2}m_{\eta^{(\prime)}}^{2}}+1\bigg{)}$
$\displaystyle\hskip
7.96674pt\times\bar{I}_{4\eta^{(\prime)}}(u)~{}+~{}\frac{2um_{\eta^{(\prime)}}^{2}(m_{c}^{2}-q^{2}-u^{2}m_{\eta^{(\prime)}}^{2})}{(m_{c}^{2}-q^{2}+u^{2}m_{\eta^{(\prime)}}^{2})^{2}}~{}~{}\bigg{(}\frac{d}{du}-\frac{6um_{\eta^{(\prime)}}^{2}}{m_{c}^{2}-q^{2}+u^{2}m_{\eta^{(\prime)}}^{2}}\bigg{)}~{}\int_{u}^{1}{d\xi}\bar{I}_{4\eta^{(\prime)}}(\xi)\bigg{]}\bigg{\\}}$
$\displaystyle\hskip
7.96674pt+\frac{m_{c}^{4}{f_{\eta^{(\prime)}}}e^{-m_{c}^{2}/M^{2}}}{4(m_{c}^{2}-q^{2}+m_{\eta^{(\prime)}}^{2})^{2}}\bigg{(}\frac{d}{du}\phi_{4;\eta}(u)\bigg{)}\bigg{|}_{u\to
1}+\frac{\alpha_{s}C_{F}e^{m_{D_{s}}^{2}/M^{2}}}{8\pi
m_{D_{s}}^{2}f_{D_{s}}}F_{1}(q^{2},m^{2},s_{0}),$

with $s(u)=m_{c}^{2}-(q^{2}-m_{\eta^{(\prime)}}^{2}u)\bar{u}/u$. The LO
invariant amplitudes include twist-2, 3, 4 contributions. The NLO QCD
corrections to the invariant amplitudes $F_{1}(q^{2},(p+q)^{2})$ include
twist-2, 3 contributions, which can be separated into the following form,

$\displaystyle F_{1}(q^{2},m^{2},s_{0})$
$\displaystyle=\frac{1}{\pi}\int_{m_{c}^{2}}^{s_{0}}dse^{-s(u)/M^{2}}{\rm
Im}F_{1}(q^{2},s)$
$\displaystyle=\frac{f_{\eta^{(\prime)}}}{\pi}\int_{m_{c}^{2}}^{s_{0}}dse^{-s(u)/M^{2}}\int_{0}^{1}{du}\Big{\\{}{\rm
Im}{T_{1}}(q^{2},s,u)\phi_{2;\eta^{(\prime)}}(u)$
$\displaystyle+\frac{\mu_{\eta^{(\prime)}}}{m_{c}}\big{[}{\rm
Im}T_{1}^{p}(q^{2},s,u)\phi_{3;\eta^{(\prime)}}^{p}(u)+{\rm
Im}T_{1}^{\sigma}(q^{2},s,u)\phi_{3;\eta^{(\prime)}}^{\sigma}(u)\big{]}\Big{\\}},$
(53)

where ${\mu_{\eta^{(\prime)}}}={m_{\eta^{(\prime)}}^{2}}/{2{m_{s}}}$. The
imaginary parts of the amplitudes $T_{1}(q^{2},s,u)$, $T_{1}^{p}(q^{2},s,u)$
and $T_{1}^{\sigma}(q^{2},s,u)$ are the hard-scattering amplitudes calculated
by the 6 diagrams for the gluon corrections. The final detailed expressions
agree with those of Refs. Duplancic:2008ix ; Duplancic:2015zna , which are not
listed here. The lower limit of integral

$\displaystyle
u_{0}=\Big{(}q^{2}-s_{0}+m_{\eta^{(\prime)}}^{2}+\sqrt{(q^{2}-s_{0}+m_{\eta^{(\prime)}}^{2})^{2}-4m_{\eta^{(\prime)}}^{2}(q^{2}-m_{c}^{2})}\Big{)}/(2m_{\eta^{(\prime)}}^{2})$
(54)

and for the final expression, we need a brief introduction to the integrals
over three-particle LCDAs, i.e.

$\displaystyle I_{3;\eta^{(\prime)}}(u)$
$\displaystyle=\frac{d}{du}\Bigg{[}\int_{0}^{u}d\alpha_{1}\int_{\Delta}^{1}dv\Phi_{3;\eta^{(\prime)}}(\alpha_{i})\Bigg{]},$
$\displaystyle I_{4;\eta^{(\prime)}}(u)$
$\displaystyle=\frac{d}{du}\Bigg{\\{}\int_{0}^{u}d\alpha_{1}\int_{\Delta}^{1}\frac{dv}{v}\Bigg{[}2\Psi_{4;\eta^{(\prime)}}(\alpha_{i})-\Phi_{4;\eta^{(\prime)}}(\alpha_{i})+2\widetilde{\Psi}_{4;\eta^{(\prime)}}(\alpha_{i})-\widetilde{\Phi}_{4;\eta^{(\prime)}}(\alpha_{i})\Bigg{]}\Bigg{\\}},$
$\displaystyle\bar{I}_{4;\eta^{(\prime)}}(u)$
$\displaystyle=\frac{d}{du}\Bigg{\\{}\int_{0}^{u}d\alpha_{1}\int_{\Delta}^{1}\frac{dv}{v}\Bigg{[}\Psi_{4;\eta^{(\prime)}}(\alpha_{i})~{}+~{}\Phi_{4;\eta^{(\prime)}}(\alpha_{i})+\widetilde{\Psi}_{4;\eta^{(\prime)}}(\alpha_{i})\,+\,\widetilde{\Phi}_{4;\eta^{(\prime)}}(\alpha_{i})\Bigg{]}\Bigg{\\}},$
(55)

where $\Delta={(u-\alpha_{1})}/{(1-\alpha_{1})}$,
$\alpha_{2}=1-\alpha_{1}-\alpha_{3}$ and $\alpha_{3}=(u-\alpha_{1})/v$. Due to
the contributions from three-particle parts are quite small, i.e. $<0.1\%$, we
can safely neglect these parts in this paper. We would like to figure out that
the decay branching fraction for the considered decay. Using the
parametrization of the transition matrix elements in terms of TFFs, in
massless lepton case, we get

$\displaystyle\frac{d\Gamma}{dq^{2}}(D_{s}^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell})$
$\displaystyle=\frac{G_{F}^{2}|V_{cs}|^{2}}{192\pi^{3}m_{D_{s}}^{3}}\Big{[}\left(m_{D_{s}}^{2}+m_{\eta^{(\prime)}}^{2}-q^{2}\right)^{2}-4m_{D_{s}}^{2}m_{\eta^{(\prime)}}^{2}\Big{]}^{3/2}|f^{\eta^{(\prime)}}_{+}(q^{2})|^{2},$
(56)

where the fermi coupling constant $G_{F}=1.166\times{10}^{-5}~{}{\rm
GeV}^{-2}$.

##  III Numerical Discussions

###  III.1 Input parameters

We adopt the following parameters to do the numerical calculation. The current
charm-quark mass is $m_{c}=1.27\pm 0.02~{}{\rm GeV}$, the masses of $D_{s}$,
$\eta$ and $\eta^{\prime}$-meson $m_{D_{s}}=1.9685~{}{\rm GeV}$,
$m_{\eta}=0.5478~{}{\rm GeV}$, $m_{\eta^{\prime}}=0.9578~{}{\rm GeV}$ and
$s$-quark mass $m_{s}=0.093~{}{\rm GeV}$. All of them are taken from the
Particle Data Group (PDG) Zyla:2020zbs . The $D_{s}$,
$\eta,\eta^{\prime}$-meson decay constants are taken as $f_{D_{s}}=0.274\pm
0.013\pm 0.007~{}{\rm GeV}$ Azizi:2010zj , $f_{\eta}=0.130\pm 0.003~{}{\rm
GeV}$ Ball:2004ye  and $f_{\eta^{\prime}}=0.157\pm 0.003~{}{\rm GeV}$
Ali:1998eb . The values of non-perturbative vacuum condensates up to
6-dimension are taken as follows Shifman:1980ui ; Colangelo:2000dp ;
Narison:2014ska ,

$\displaystyle\langle\alpha_{s}G^{2}\rangle=0.038\pm 0.011~{}{\rm GeV}^{4},$
$\displaystyle\langle g_{s}^{3}fG^{3}\rangle=0.045\pm 0.007~{}{\rm GeV}^{6},$
$\displaystyle\langle g_{s}\bar{q}q\rangle^{2}=(2.082_{-0.697}^{+0.734})\times
10^{-3}~{}{\rm GeV}^{6},$ $\displaystyle\langle
g_{s}^{2}\bar{q}q\rangle^{2}=(7.420_{-2.483}^{+2.614})\times 10^{-3}~{}{\rm
GeV}^{6},$ $\displaystyle\langle
q\bar{q}\rangle=(-2.417_{-0.114}^{+0.227})\times 10^{-2}~{}{\rm GeV}^{3},$
$\displaystyle\langle\bar{s}s\rangle=(-1.789_{-0.084}^{+0.168})\times
10^{-2}~{}{\rm GeV}^{3},$ $\displaystyle\langle
g_{s}s\bar{s}\rangle^{2}=(1.541_{-0.516}^{+0.543})\times{10^{-3}}~{}{\rm
GeV}^{6},$ $\displaystyle\sum{\langle
g_{s}^{2}\bar{q}q\rangle^{2}}=(1.891_{-0.633}^{+0.665})\times{10^{-2}}~{}{\rm
GeV}^{6}.$ (57)

The quark-gluon mixture condensate $\langle g_{s}\bar{q}\sigma
TGq\rangle=m_{0}^{2}\langle\bar{q}q\rangle$ with $m_{0}^{2}=0.80\pm 0.02{\rm
GeV}^{2}$, which leads to

$\displaystyle\langle g_{s}\bar{s}\sigma
TGs\rangle=(-1.431_{-0.076}^{+0.139})\times{10^{-2}}~{}{\rm GeV}^{5}.$ (58)

Here, the ratio $\kappa=\langle\bar{s}s\rangle/\langle\bar{q}q\rangle=0.74\pm
0.03$ has been used. Meanwhile, the typical scale in this paper is $\mu_{\rm
IR}=(m_{D_{s}}^{2}-m_{c}^{2})^{1/2}\approx 1.5~{}{\rm GeV}$. So the
renormalization group equations (RGEs) should be used for running the quark
masses and each vacuum condensates appearing in the BFTSR from the initial
scale $\mu_{0}=1~{}{\rm GeV}$ to the typical scale $\mu_{\rm IR}$. The RGE can
be found in Refs. Yang:1993bp ; Hwang:1994vp ; Lu:2006fr , which are not
listed here.

Table 1: The criteria for determining Borel windows, and the resultant Borel windows and the corresponding values of the $\eta$ and $\eta^{\prime}$-meson leading twist LCDA moments $\langle\xi_{2;\eta^{(\prime)}}^{n}\rangle|_{\mu}$. “Con.” represents the continuum contribution and “Six.” represents the dimension-six condensates’ contribution. | $n$ |  $0$ |  $2$ |  $4$ | $6$
---|---|---|---|---|---
| Con. | $<20\%$ | $<25\%$ | $<30\%$ | $<35\%$
| Six. | $<5\%$ | $<5\%$ | $<5\%$ | $<5\%$
$\eta$-meson | ${M^{2}}$ | $[0.535,1.188]$ | $[1.026,1.402]$ | $[1.368,1.759]$ | $[1.677,2.194]$
| $\langle{\xi^{n}_{2;\eta}}\rangle|_{\mu}$ | $[0.952,1.168]$ | $[0.231,0.230]$ | $[0.110,0.102]$ | $[0.067,0.059]$
| $n$ |  $0$ |  $2$ |  $4$ | $6$
| Con. | $<30\%$ | $<35\%$ | $<40\%$ | $<45\%$
| Six. | $<5\%$ | $<5\%$ | $<5\%$ | $<5\%$
$\eta^{\prime}$-meson | ${M^{2}}$ | $[1.049,1.137]$ | $[1.026,1.233]$ | $[1.368,1.627]$ | $[1.677,2.082]$
| $\langle{\xi^{n}_{2;\eta^{\prime}}}\rangle|_{\mu}$ | $[1.076,1.061]$ | $[0.221,0.201]$ | $[0.099,0.086]$ | $[0.059,0.047]$

###  III.2 Determination for the Gegenbauer moments of $\eta^{(\prime)}$-meson
twist-2 LCDA

One of the significant parameters in BFTSR is the continuum threshold
$s_{\eta^{(\prime)}}$ for the moments of $\eta^{(\prime)}$-meson twist-2
LCDAs. Following our previous works, we can determine $s_{\eta}=1.3\pm
0.1~{}{\rm GeV}$ and $s_{\eta^{\prime}}=0.8\pm 0.1~{}{\rm GeV}$ by setting the
0th-order of Gegenbauer moment into 1, i.e.
$\langle\xi_{2;\eta}^{0}\rangle|_{\mu}=\langle\xi_{2;\eta^{\prime}}^{0}\rangle|_{\mu}=1$.
Meanwhile, in order to determine the allowable range of the Borel parameter
$M^{2}$ (i.e. the Borel Window), we adopt the following three criteria

  * •

The continuum contributions are less than $45\%$ of the total results;

  * •

The contributions from the dimension-six condensates do not exceed $5\%$;

  * •

We require the variations of
$\langle\xi_{2;\eta^{(\prime)}}^{n}\rangle|_{\mu}$ within the Borel window to
be less than $10\%$.

Figure 1: Contributions from the continuum state and dimension-six condensates
for the $\eta^{(\prime)}$-meson leading-twist LCDA moments
$\langle{\xi^{n}_{2;\eta^{(\prime)}}}\rangle|_{\mu}$ versus the Borel
parameter $M^{2}$, where all input parameters are set to be their central
values.

For the first four Gegenbauer moments
$\langle\xi_{2;\eta^{(\prime)}}^{n}\rangle|_{\mu}$ with $n=(0,2,4,6)$, we list
the allowable Borel region and their corresponding
$\langle\xi_{2;\eta^{(\prime)}}^{n}\rangle|_{\mu}$ in Table 1. When
$n=(0,2,4,6)$, we have set the continuum contributions to be less than $20\%$,
$25\%$, $30\%$, $35\%$ for $\langle\xi_{2;\eta}^{n}\rangle|_{\mu}$ and $30\%$,
$35\%$, $40\%$, $45\%$ for $\langle\xi_{2;\eta^{\prime}}^{n}\rangle|_{\mu}$,
respectively, and the dimension-six condensates’ contributions to be less than
$5\%$ for all the order of $\langle\xi_{2;\eta^{(\prime)}}^{n}\rangle|_{\mu}$.
Then, the determined continuum’s and dimension-six condensates’ contribution
for $\langle\xi_{2;\eta}^{n}\rangle|_{\mu}$ and
$\langle\xi_{2;\eta^{\prime}}^{n}\rangle|_{\mu}$ with $n=(2,4,6)$ are shown in
the left and right panel of Fig. 1, respectively. In which, the shaded region
stand for the Borel windows.

In order to provide a deeper insight into the flatness of the LCDA moments
versus the Borel parameter $M^{2}$, we present the first three curves for the
moments of $\eta^{(\prime)}$-meson twist-2 LCDA at the initial scale, i.e.
$\langle\xi^{n}_{2;\eta^{(\prime)}}\rangle|_{\mu}$ with $n=(2,4,6)$ in Figure
2. The determined Borel window are $M^{2}\in[1.0,2.5]~{}{\rm GeV}^{2}$ for
$\eta$-meson and $M^{2}\in[1.0,2.4]~{}{\rm GeV}^{2}$ for
$\eta^{\prime}$-meson. It is noted that within this range, the moments
$\langle\xi^{2}_{2;\eta^{(\prime)}}\rangle|_{\mu}$,
$\langle\xi^{4}_{2;\eta^{(\prime)}}\rangle|_{\mu}$ and
$\langle\xi^{6}_{2;\eta^{(\prime)}}\rangle|_{\mu}$ are almost flat, which vary
less than $10\%$ for the total results in the Borel window.

By taking the squared average of all the uncertainty sources into
consideration and making use of the relations between the Gegenbauer moments
$a^{n}_{2;\eta^{(\prime)}}(\mu)$ and the LCDA moments
$\langle\xi^{n}_{2;\eta^{(\prime)}}\rangle|_{\mu}$, i.e. Eq. (43), we obtain
the first three $a^{n}_{2;\eta^{(\prime)}}(\mu_{0})$ and
$\langle\xi^{n}_{2;\eta^{(\prime)}}\rangle|_{\mu_{0}}$ with $n=(2,4,6)$ for
the leading-twist $\eta^{(\prime)}$-meson LCDA
$\phi_{2;\eta^{(\prime)}}(u,\mu_{0})$ predicted from BFTSR in Table 2. The
factorization scale is taken as the initial scale ${\mu_{0}}=1~{}{\rm GeV}$.
As a comparison, we also list the LCSR given in year 2013 Offen:2013nma , the
CLEO fit Gronberg:1997fj , the BABAR fit BABAR:2011ad , the ones given by
Kroll Kroll:2013iwa  and Ball Ball:2004ye , respectively. For the Ball’s
results, it is calculated by using the approximation
$a_{2;\eta}^{n}(\mu)=a_{2;K}^{n}(\mu)=a_{2;\pi}^{n}(\mu)$. Our results for the
second and fourth order $\eta$-meson LCDA’s moments, e.g.
$\langle\xi_{2;\eta}^{2}\rangle|_{\mu_{0}}$ and
$\langle\xi_{2;\eta}^{4}\rangle|_{\mu_{0}}$ agree with the Ball’s predictions
within errors. But there still exist discrepancy for $a_{2;\eta}^{4}(\mu_{0})$
with Ball’s prediction. The main reason lies in the fourth order of equations
between Gegenbauer and LCDA moments, i.e. Eq. (43) have large coefficients
which will enlarge the small discrepancy of LCDA moments, even appears the
opposite sign. At present, there are few studies on $\eta^{\prime}$-meson’s
twist-2 LCDA.

Figure 2: The first three moments $\langle\xi^{n}_{2;\eta^{(\prime)}}\rangle|_{\mu}$ with $(n=2,4,6)$ versus the Borel parameter $M^{2}$. The darker shaded bands indicate the Borel windows for $\langle\xi^{n}_{2;\eta^{(\prime)}}\rangle|_{\mu}$, respectively. Table 2: First three Gegenbaner and LCDA moments $a^{n}_{2;\eta^{(\prime)}}(\mu_{0})$ and $\langle\xi^{n}_{2;\eta^{(\prime)}}\rangle|_{\mu_{0}}$ with $n=(2,4,6)$ for the leading-twist $\eta^{(\prime)}$-meson LCDA $\phi_{2;\eta^{(\prime)}}(u,\mu_{0})$, where the errors are squared averages of those from all the input parameters. Other theoretical predictions are also given as a comparison. |  |  $\eta$-meson |  |  |  |  $\eta^{\prime}$-meson |
---|---|---|---|---|---|---|---
| $\langle\xi^{2}_{2;\eta}\rangle|_{\mu_{0}}$ | $\langle\xi^{4}_{2;\eta}\rangle|_{\mu_{0}}$ | $\langle\xi^{6}_{2;\eta}\rangle|_{\mu_{0}}$ |  | $\langle\xi^{2}_{2;\eta^{\prime}}\rangle|_{\mu_{0}}$ | $\langle\xi^{4}_{2;\eta^{\prime}}\rangle|_{\mu_{0}}$ | $\langle\xi^{6}_{2;\eta^{\prime}}\rangle|_{\mu_{0}}$
BFTSR (This work) | $0.231_{-0.013}^{+0.010}$ | $0.109_{-0.007}^{+0.007}$ | $0.066_{-0.006}^{+0.006}$ |  | $0.211_{-0.017}^{+0.015}$ | $0.093_{-0.009}^{+0.009}$ | $0.054_{-0.008}^{+0.008}$
CLEO fit Gronberg:1997fj  | $0.176\pm 0.010$ | - | - |  | - | - | -
BABAR fit BABAR:2011ad  | $0.183\pm 0.007$ | - | - |  | - | - | -
P. Kroll Kroll:2013iwa  | $0.183\pm 0.007$ | - | - |  | - | - | -
SR fit Offen:2013nma  | $0.286\pm 0.051$ | - | - |  | - | - | -
P. Ball Ball:2004ye  | 0.239 | 0.110 | - |  | - | - | -
|  |  $\eta$-meson |  |  |  |  $\eta^{\prime}$-meson |
| $a^{2}_{2;\eta}(\mu_{0})$ | $a^{4}_{2;\eta}(\mu_{0})$ | $a^{6}_{2;\eta}(\mu_{0})$ |  | $a^{2}_{2;\eta^{\prime}}(\mu_{0})$ | $a^{4}_{2;\eta^{\prime}}(\mu_{0})$ | $a^{6}_{2;\eta^{\prime}}(\mu_{0})$
BFTSR (This work) | $0.090_{-0.037}^{+0.031}$ | $0.025_{-0.010}^{+0.003}$ | $0.033_{-0.058}^{+0.055}$ |  | $0.033_{-0.050}^{+0.042}$ | $-0.002_{-0.016}^{+0.007}$ | $0.043_{-0.072}^{+0.067}$
CLEO fit Gronberg:1997fj  | $-0.07\pm 0.03$ | - | - |  | - | - | -
BABAR fit BABAR:2011ad  | $-0.05\pm 0.02$ | - | - |  | - | - | -
P. Kroll Kroll:2013iwa  | $-0.05\pm 0.02$ | - | - |  | - | - | -
SR fit Offen:2013nma  | $0.25\pm 0.15$ | - | - |  | - | - | -
P. Ball Ball:2004ye  | $0.115$ | $-0.015$ | - |  | - | - | -

After considering the Gegenbauer moments $a_{2;\eta^{(\prime)}}^{n}(\mu)$ up
to 6th-order into the conformal expansion of the Gegenbauer polynomial at
initial scale, i.e. Eq. (10), we present the $\eta$ and $\eta^{\prime}$-meson
twist-2 LCDAs in Figure 3(a) and Figure 3(b) separately. For
$\phi_{2;\eta}(u,\mu_{0})$, we present the asymptotic form, the CLEO
Gronberg:1997fj , the BABAR BABAR:2011ad , Kroll’s prediction Kroll:2013iwa ,
the SR fit Offen:2013nma  and the Ball’s prediction Ball:2004ye  as a
comparison. Figure 3(a) shows that the LCSR 2013 and Ball’s results prefer a
double-peaked behavior. The reason lies in that they adopt the $\pi,K$-meson’s
LCDAs as those of $\eta$-meson LCDAs. Conversely, the CLEO Gronberg:1997fj ,
the BABAR BABAR:2011ad  and the fitting results by Kroll indicate a single-
peaked behavior. Our prediction tends to a double-peaked behavior. For
$\phi_{2;\eta^{\prime}}(u,\mu_{0})$, we only exhibit the asymptotic form due
to there are less results from references, which is shown in Figure 3(b).
Furthermore, in order to have a look at the evolution
$\phi_{2;\eta^{(\prime)}}(u,\mu_{0})$ with $n=(2,4,6)$, we present the
different curves in Figure 4. If we take $n=(2,4)$, the behavior of LCDAs
shall be closer to the asymptotic form. Furthermore, a small shake is observed
when taking the 6th-order LCDA moment into consideration.

Figure 3: The $\eta^{(\prime)}$-meson leading-twist LCDA
$\phi_{2;\eta^{(\prime)}}(u,{\mu}_{0})$ predicted from the BFTSR. We make a
comparison with the asymptotic form, the CLEO Gronberg:1997fj , SR fit
Offen:2013nma , the BABAR BABAR:2011ad , and the predictions of Kroll
Kroll:2013iwa  and Ball Ball:2004ye .

Figure 4: The curves of $\eta$ and $\eta^{\prime}$-meson twist-2 LCDA with
$n=(2,4,6)$ respectively.

Other two-particle Fock state twist-3 and twist-4 LCDAs
$\phi_{3;\eta^{(\prime)}}^{p}(u)$, $\phi_{3;\eta^{(\prime)}}^{\sigma}(u)$,
${\psi_{4;\eta^{(\prime)}}}(u)$ and ${\phi_{4;\eta^{(\prime)}}}(u)$ are
defined as follows

$\displaystyle\phi_{3;\eta^{(\prime)}}^{p}(u)$
$\displaystyle=1+\Big{(}30\eta_{3}^{\eta^{(\prime)}}-\frac{5}{2}\rho_{\eta^{(\prime)}}^{2}\Big{)}C_{2}^{1/2}(\xi)+\Big{(}-3\eta_{3}^{\eta^{(\prime)}}\omega_{3}^{\eta^{(\prime)}}-\frac{{27}}{{20}}\rho_{\eta^{(\prime)}}^{2}-\frac{81}{10}\rho_{\eta^{(\prime)}}^{2}a^{2}_{2;\eta^{(\prime)}}\Big{)}$
$\displaystyle\times C_{4}^{1/2}(\xi),$ (59)
$\displaystyle\phi_{3;\eta^{(\prime)}}^{\sigma}(u)$
$\displaystyle=6u\bar{u}\Big{(}1+5\eta_{3}^{\eta^{(\prime)}}-\frac{1}{2}\eta_{3}^{\eta^{(\prime)}}\omega_{3}^{\eta^{(\prime)}}-\frac{7}{{20}}\rho_{\eta^{(\prime)}}^{2}-\frac{3}{5}\rho_{\eta^{(\prime)}}^{2}a^{2}_{2;\eta^{(\prime)}}\Big{)}C_{2}^{3/2}(\xi),$
(60) $\displaystyle\psi_{4;\eta^{(\prime)}}(u)$
$\displaystyle=\frac{5}{2}\varepsilon^{2}u^{2}\bar{u}^{2}+\frac{1}{2}\varepsilon\delta^{2}\bigg{[}u\bar{u}(2+13u\bar{u})+10u^{3}\ln
u\ln\bar{u}(2-3u+\frac{6}{5}u^{2})+10\,\bar{u}^{3}\,(2$
$\displaystyle-3\bar{u}+\frac{6}{5}{\bar{u}^{2}})\bigg{]},$ (61)
$\displaystyle\phi_{4;\eta^{(\prime)}}(u)$
$\displaystyle=\frac{{10}}{3}{\delta^{2}}u\bar{u}(u-\bar{u}),$ (62)

where, the values for the twist-3, 4 LCDAs parameters are taken from Refs.
Ball:2006wn ; Huang:2001xb . In order to run the hadronic parameters of the
$\eta^{(\prime)}$-meson twist-2, 3, 4 LCDAs from the initial factorization
scale to any other scale, especially for typical scale $\mu_{\rm IR}$, the
renormalization group equation should be used, which has the form

$\displaystyle c_{i}(\mu_{\rm IR})={\cal
L}^{\gamma_{c_{i}}/\beta_{0}}c_{i}(\mu_{0}),$ (63)

where ${\cal L}=\alpha_{s}(\mu_{\rm IR})/\alpha_{s}(\mu_{0})$,
$\beta_{0}=11-2/3n_{f}$, and the one-loop anomalous dimensions
$\gamma_{c_{i}}$ can be seen in our previous work Fu:2020uzy . Taking the
hadronic parameters at initial scale $\mu_{0}$ and using the renormalization
function (63), one can achieve the corresponding values at the typical scale
$\mu_{\rm IR}$.

Table 3: The TFFs $f_{+}^{\eta^{(\prime)}}(0)$ at the large recoil point $q^{2}=0$. As a comparison, we also present the predictions from various experimental and theoretical groups. | References | $f_{+}^{\eta}(0)$ | $f_{+}^{\eta^{\prime}}(0)$
---|---|---|---
| BESIII Ablikim:2019rjz  | $0.4576(70)$ | $0.490(51)$
Experimental results | LQCD-I Bali:2014pva  | $0.542(13)$ | $0.404(25)$
| LQCD-II Bali:2014pva  | $0.564(11)$ | $0.437(18)$
| This work (LCSR) | $0.476_{-0.036}^{+0.040}$ | $0.544_{-0.042}^{+0.046}$
| LFQM Verma:2011yw  | $0.76$ | -
| CQM Melikhov:2000yu  | $0.78$ | $0.78$
| CCQM Soni:2018adu  | $0.78(12)$ | $0.73(11)$
Theoretical predictions | CCQM Ivanov:2019nqd  | $0.49(7)$ | $0.59(9)$
| LCSR 2013 Offen:2013nma  | $0.432(33)$ | $0.520(80)$
| LCSR 2015 Duplancic:2015zna  | $0.495_{-0.030}^{+0.029}$ | $0.558_{-0.045}^{+0.047}$
| QCD SR Colangelo:2001cv  | $0.50(4)$ | -
Figure 5: The TFFs $f_{+}^{\eta^{(\prime)}}(0)$ versus the Borel parameters
$M^{2}$, where the shaded band is induced by the variations of squared average
of all input parameters.

###  III.3 TFFs and series expansion

In order to determine the continuum threshold $s_{0}$ for the
$D_{s}\to\eta^{(\prime)}$ TFFs within LCSR approach, i.e. Eq. (LABEL:Eq:fp),
one can follow the four criteria

  * •

The continuum contributions are less than $30\%$ of the total results;

  * •

The contributions from the twist-4 LCDAs do not exceed $5\%$;

  * •

We require the variations of the TFF within the Borel window be less than
$10\%$;

  * •

The continuum threshold $s_{0}$ should be closer to the squared mass of the
first excited state of $D_{s}$-meson.

Based on the fourth term of the criteria, we take $s_{0}$ to be close to the
squared mass of the excited state of $D_{s}$-meson $D_{s1}(2460)$, i.e.
$s_{0}=6.1(3)~{}{\rm GeV^{2}}$. Furthermore, the Borel parameter is taken as
$M^{2}=25(1)~{}{\rm GeV^{2}}$. Furthermore, we can obtain the sum rule for
$m_{D_{s}}$ by differentiating the form factors
$f_{+}^{\eta^{(\prime)}}(q^{2})f_{D_{s}}$ with respect to $-1/M^{2}$
Fu:2014pba . The resultant $m_{D_{s}}^{\rm LCSR}=1.9653~{}{\rm GeV}$ agrees
with the measured value $m_{D_{s}}^{\rm PDG}=1.9685~{}{\rm GeV}$. In order to
show the degree of stability of the TFFs versus the Borel parameter, we
present the curve of TFFs in Figure 5, in which the shaded region shows the
errors from all input parameters. The solid line with blue shaded band
represents $f_{+}^{\eta}(0)$ and the dashed line with green band represents
$f_{+}^{\eta^{\prime}}(0)$. This figure shows that
$f_{+}^{\eta^{(\prime)}}(0)$ changes less than $0.5\%$ within the range of
$M^{2}\in[23,27]~{}{\rm GeV^{2}}$, which satisfies the third term of the
criteria. More definitely, we put the errors caused by different input
parameters in the following,

$\displaystyle f_{+}^{\eta}(0)$
$\displaystyle=0.476+(_{-0.012}^{+0.011})_{s_{0}}+(_{-0.001}^{+0.001})_{M^{2}}+(_{-0.032}^{+0.036})_{m_{c},f_{D_{s}}}+(_{-0.011}^{+0.011})_{f_{\eta}}+(_{-0.003}^{+0.004})_{a_{2;\eta}^{2}}$
$\displaystyle+(_{-0.000}^{+0.000})_{a_{2;\eta}^{4}}+(_{-0.000}^{+0.000})_{a_{2;\eta}^{6}}$
$\displaystyle=0.476_{-0.036}^{+0.040}$ (64) $\displaystyle
f_{+}^{\eta^{\prime}}(0)$
$\displaystyle=0.544+(_{-0.016}^{+0.015})_{s_{0}}+(_{-0.001}^{+0.001})_{M^{2}}+(_{-0.037}^{+0.042})_{m_{c},f_{D_{s}}}+(_{-0.010}^{+0.010})_{f_{\eta^{\prime}}}+(_{-0.005}^{+0.006})_{a_{2;\eta^{\prime}}^{2}}$
$\displaystyle+(_{-0.001}^{+0.000})_{a_{2;\eta^{\prime}}^{4}}+(_{-0.000}^{+0.000})_{a_{2;\eta^{\prime}}^{6}}$
$\displaystyle=0.544_{-0.042}^{+0.046}$ (65)

According to Eq. (10), when the Gegenbauer moments are taken up to 2nd, 4th
and 6th order levels, the central values of TFFs $f_{+}^{\eta^{(\prime)}}(0)$
are

$\displaystyle f_{+}^{\eta}(0)|_{\rm 2nd}=0.4745,$ $\displaystyle
f_{+}^{\eta}(0)|_{\rm 4th}=0.4755,$ $\displaystyle f_{+}^{\eta}(0)|_{\rm
6th}=0.4763,$ $\displaystyle f_{+}^{\eta^{\prime}}(0)|_{\rm 2nd}=0.5435,$
$\displaystyle f_{+}^{\eta^{\prime}}(0)|_{\rm 4th}=0.5434,$ $\displaystyle
f_{+}^{\eta^{\prime}}(0)|_{\rm 6th}=0.5436.$

In comparing with the 2nd-order Gegenbauer moment’s contribution, the 4th,
6th-order contributions shall be changed by about 0.211%, 0.379% for
$f_{+}^{\eta}(0)$ and $-0.018\%$, 0.018% for $f_{+}^{\eta^{\prime}}(0)$,
respectively. These ratios are really small, indicating the Gegenbauer series
has good convergence over the moment expansion. Then, we list the TFFs for
$D_{s}\to\eta^{(\prime)}$ at large recoil point, i.e.
$f_{+}^{\eta^{(\prime)}}(0)$, in Table 3, in which the uncertainties are from
the squared average of all the mentioned error sources. As a comparison, we
also present other theoretical and experimental predictions, such as the
BESIII Ablikim:2019rjz , the LQCD Bali:2014pva , the LFQM Verma:2011yw , the
CQM Melikhov:2000yu , the CCQM Soni:2018adu ; Ivanov:2019nqd , the QCD SR
Colangelo:2001cv , the LCSR at 2013 and 2015 Offen:2013nma ; Duplancic:2015zna
, respectively. Our results agree with the BESIII, the CCQM, the LCSR, the
LQCD within errors, but are lack of agreement with the Lattice QCD results.

The physically allowable ranges for the TFFs are $m_{\ell}^{2}\leq
q^{2}\leq(m_{D_{s}}-m_{\eta})^{2}\approx 2~{}{\rm GeV^{2}}$ and
$m_{\ell}^{2}\leq q^{2}\leq(m_{D_{s}}-m_{\eta^{\prime}})^{2}\approx 1~{}{\rm
GeV^{2}}$. Theoretically, the LCSRs approach for $D_{s}\to\eta^{(\prime)}$
TFFs are applicable in low and intermediate $q^{2}$-regions, i.e.
$q^{2}\in[0,1.2]~{}{\rm GeV^{2}}$ of $\eta$-meson, $q^{2}\in[0,0.6]~{}{\rm
GeV^{2}}$ of $\eta^{\prime}$-meson. One can extrapolate it to whole
$q^{2}$-regions via a rapidly $z(q^{2},t)$ converging series expansion (SE)
Bharucha:2010im :

$\displaystyle
f^{\eta^{(\prime)}}_{+}(q^{2})=\frac{1}{B(q^{2})\phi_{+}^{\eta^{(\prime)}}(q^{2})}\sum\limits_{k}^{K-1}\alpha_{k}z^{k}(q^{2},t_{0}).$
(66)

Table 4: The fitting parameters $z(q^{2},t_{0}^{\eta^{(\prime)}})$ of SSE for TFFs $f_{+}^{{\eta^{(\prime)}}}(q^{2})$. | $q^{2}({\rm GeV^{2}})$ | 0 | 0.4 | 0.8 | 1.2 | 1.6 | 2.0
---|---|---|---|---|---|---|---
$\eta$-meson | $z(q^{2},t_{0}^{\eta})$ | $0.048$ | $0.032$ | $0.014$ | $-0.005$ | $-0.025$ | $-0.047$
| $q^{2}({\rm GeV^{2}})$ | 0 | 0.2 | 0.4 | 0.6 | 0.8 | 1.0
$\eta^{\prime}$-meson | $z(q^{2},t_{0}^{\eta^{\prime}})$ | $0.016$ | $0.004$ | $-0.009$ | $-0.022$ | $-0.036$ | $-0.051$

Normally, the parameter $K$ stands for the order of expansion. The
$z(q^{2},t)$ is the function:

$\displaystyle
z(q^{2},t)=\frac{\sqrt{t_{+}-q^{2}}-\sqrt{t_{+}-t}}{\sqrt{t_{+}-q^{2}}+\sqrt{t_{+}-t}},$
(67)

where $t=(t_{0}^{\eta^{(\prime)}},t_{-}^{\eta^{(\prime)}},m_{D_{s}})$ with
$t_{\pm}^{\eta^{(\prime)}}=(m_{D_{s}}\pm m_{\eta^{(\prime)}})^{2}$. Here, the
$0\leq t_{0}^{\eta^{(\prime)}}\leq t_{-}^{\eta^{(\prime)}}$ is a free
parameter which can be optimised to reduce the maximum value of
$|z(q^{2},t_{0}^{\eta^{(\prime)}})|$ in the physical TFFs range,
$t_{0}^{\eta^{(\prime)}}|_{\rm
opt.}=t_{+}^{\eta^{(\prime)}}(1-\sqrt{1-t_{-}^{\eta^{(\prime)}}/t_{+}^{\eta^{(\prime)}}})$.
Here, we list the value of $z(q^{2},t_{0}^{\eta^{(\prime)}})$ with different
$q^{2}$ cases in Table 4. The function $\phi_{+}^{\eta^{(\prime)}}(q^{2})$ can
be expressed as

$\displaystyle\phi_{+}^{\eta^{(\prime)}}(q^{2})$ $\displaystyle=$
$\displaystyle\sqrt{\dfrac{\varsigma}{48\pi\chi_{+}^{\eta^{(\prime)}}(n)}}\dfrac{q^{2}-t_{+}}{(t_{+}-t_{0})^{1/4}}\bigg{(}\dfrac{z(q^{2},0)}{-q^{2}}\bigg{)}^{(3+n)/2}\bigg{(}\dfrac{z(q^{2},t_{0})}{t_{0}-q^{2}}\bigg{)}^{-1/2}$
(68) $\displaystyle\times$
$\displaystyle\bigg{(}\dfrac{z(q^{2},t_{-})}{t_{-}-q^{2}}\bigg{)}^{-3/4}\bigg{|}_{n=2},$

Figure 6: The SSE for the TFFs $f^{\eta}_{+}(q^{2})$ and
$f^{\eta^{\prime}}_{+}(q^{2})$ up to $K=3,4,5$ order.

where $\varsigma=1$ is an isospin-degeneracy factor for a given channel and
$B(q^{2})=\prod\limits_{i}z(q^{2},m_{D_{s}^{i}}^{2})$ stands for the Blaschke
factor. The $m_{D_{s}^{i}}$ stands for the mass of each resonance state with
$J^{P}=0^{-},0^{+},1^{-}...$, respectively, which can be found in PDG
Zyla:2020zbs . For the coefficients $\chi_{+}^{\eta^{(\prime)}}(n)$, it can be
calculated by using QCD sum rules including perturbative LO and NLO results as
well as the condensate contributions, which can be expressed as
Bharucha:2010im

$\displaystyle\chi_{+}^{\eta^{(\prime)}}=\frac{3}{32\pi^{2}m_{c}^{2}}\bigg{[}1+\frac{\alpha_{s}(m_{c})C_{F}}{4\pi}\bigg{(}\frac{25}{6}+\frac{2\pi^{2}}{3}\bigg{)}\bigg{]}-\frac{\langle
q\bar{q}\rangle}{m_{c}^{5}}-\frac{\langle\alpha_{s}G^{2}\rangle}{12\pi
m_{c}^{6}}-\frac{\langle\bar{q}Gq\rangle}{m_{c}^{7}}.$ (69)

Furthermore, the coefficients $\alpha_{k}$ should satisfy the basic unitarity
constraint,

$\displaystyle\sum\limits_{k}^{K-1}\alpha_{k}^{2}<1$ (70)

The simplified version of the series expansion (SSE) method is to replace the
Blaschke factor $B(q^{2})$ by a simple pole $P(q^{2})$ to account for low-
lying resonances, i.e. Bourrely:2008za

$\displaystyle
f_{+}^{\eta^{(\prime)}}(q^{2})=\frac{1}{P(t)}\sum\limits_{k}^{K-1}\beta_{k}z^{k}(q^{2},t_{0}).$
(71)

As for the SSE parameterization, imposing the unitarity bound and by comparing
the SE and SSE parameterizations

$\displaystyle\alpha_{i}=\sum\limits_{k=0}^{\min[K-1,i]}\zeta_{i-k}\beta_{k},~{}~{}~{}~{}~{}~{}~{}0\leq
i\leq K-1,$ (72)

we obtain the unitarity bound of the SSE

$\displaystyle\sum\limits_{j,k=0}^{\min[K-1]}C_{jk}\beta_{j}\beta_{k}\leq 1,$
(73)

with the positive defined matrix

$\displaystyle C_{jk}=\sum\limits_{i=0}^{K-1-\max[j,k]}\zeta_{i}\zeta_{j-k}.$
(74)

Table 5: The fitting parameters for the central TFFs $f_{+}^{{\eta^{(\prime)}}(C)}(q^{2})$, the upper TFFs $f_{+}^{{\eta^{(\prime)}}(U)}(q^{2})$ and the lower TFFs $f_{+}^{{\eta^{(\prime)}}(L)}(q^{2})$. | $f_{+}^{\eta(C)}(q^{2})$ | $f_{+}^{\eta(U)}(q^{2})$ | $f_{+}^{\eta(L)}(q^{2})$ | $f_{+}^{\eta^{\prime}(C)}(q^{2})$ | $f_{+}^{\eta^{\prime}(U)}(q^{2})$ | $f_{+}^{\eta^{\prime}(L)}(q^{2})$
---|---|---|---|---|---|---
$\alpha_{0}$ | $-0.0003$ | $-0.0003$ | $-0.0003$ | $-0.00001$ | $-0.00001$ | $-0.00001$
$\alpha_{1}$ | $-0.0015$ | $-0.0017$ | $-0.0015$ | $-0.00017$ | $-0.00018$ | $-0.00016$
$\alpha_{2}$ | $-0.0061$ | $-0.0074$ | $-0.0052$ | $-0.00134$ | $-0.00143$ | $-0.00127$
$\displaystyle\sum{\alpha_{k}^{2}}$ | $3.9\times 10^{-5}$ | $5.8\times 10^{-5}$ | $2.9\times 10^{-5}$ | $1.8\times 10^{-6}$ | $2.4\times 10^{-6}$ | $1.6\times 10^{-6}$
$\beta_{0}$ | $0.512$ | $0.555$ | $0.476$ | $0.563$ | $0.611$ | $0.519$
$\beta_{1}$ | $-1.450$ | $-1.211$ | $-0.976$ | $-1.653$ | $-1.801$ | $-1.534$
$\beta_{2}$ | $17.26$ | $12.24$ | $8.73$ | $28.95$ | $30.61$ | $27.98$
$\displaystyle\sum C_{i,j}\beta_{i}\beta_{j}$ | $0.020$ | $0.019$ | $0.008$ | $0.005$ | $0.005$ | $0.004$
$\Delta$ | $0.075\%$ | $0.129\%$ | $0.099\%$ | $0.005\%$ | $0.004\%$ | $0.007\%$

Theoretically, the order $K$ can be taken up to infinite order. It has been
proven that the higher order expansion shall give the same result with very
small errors. We take $K=3,4,5$ as explicit examples, which are shown in
Figure 6. It is found that three SSE curves with $K=3,4,5$ for the TFFs
$f_{+}^{\eta^{(\prime)}}(q^{2})$ are almost coincide with each other. So, we
will take $K=3$ to do our expansion in the following calculations, which also
agrees with the choices of most of the theoretical groups. At the same time,
the $\alpha_{k}$ of SE and $\beta_{k}$ of SSE should also satisfy the
condition $\Delta<1\%$. From which, the parameter $\Delta$ is used to measure
the quality of extrapolation, which is defined as

$\displaystyle\Delta=\dfrac{\sum_{t}|F_{i}(t)-F_{i}^{\rm
fit}(t)|}{\sum_{t}|F_{i}(t)|}\times 100,$ (75)

where $t\in[0,1/100,\cdots,100/100]\times 1.2~{}{\rm GeV}$ for $\eta$-meson,
$t\in[0,1/100,\cdots,100/100]\times 0.6~{}{\rm GeV}$ for
$\eta^{\prime}$-meson. Numerical values of $\alpha_{k}$, $\beta_{k}$ with
$k=(0,1,2)$ for central values, upper and lower limit of
$f_{+}^{\eta^{(\prime)}}(q^{2})$ are listed in Table 5. Here, the
$\sum\alpha_{k}^{2}$ for SE unitarity bound and $\displaystyle\sum
C_{i,j}\beta_{i}\beta_{j}$ for SSE are less than 1. The quality of
extrapolation $\Delta$ are less than 0.13%. Since there have good unitary
bound and a small $\Delta$ value, the SSE results are in high agreement with
our LCSR results.

Figure 7: The TFFs $f^{\eta^{(\prime)}}_{+}(q^{2})$ together with its
uncertainties. Results of the LCSR 2013 Offen:2013nma , the LCSR 2015
Duplancic:2015zna , the CCQM Ivanov:2019nqd  and the BESIII Ablikim:2019rjz
are also given as the comparison.

The extrapolated TFFs in the whole $q^{2}$-region is shown in Figure 7, other
theoretical and experimental results, such as those of the BESIII
Ablikim:2019rjz , the LCSR 2013 Offen:2013nma , the LCSR 2015
Duplancic:2015zna  and the CCQM Ivanov:2019nqd  are present as a comparison.
Two sets of BESIII are from the two different $\eta^{(\prime)}$ decay
channels. For $\eta$-meson, the blue triangle stands for $\eta\to\gamma\gamma$
channel and the red diamond stands for $\eta\to\pi^{0}\pi^{+}\pi^{-}$ channel.
For $\eta^{\prime}$-meson the blue triangle stands for
$\eta^{\prime}\to\gamma\rho^{0}$ channel and red diamond stands for
$\eta^{\prime}\to\eta_{\gamma\gamma}\pi^{+}\pi^{-}$ channel. To compare with
other theoretical and experimental groups, our results have the following
characteristics,

  * •

Comparing with $f^{\eta}_{+}(q^{2})$, the $f^{\eta^{\prime}}_{+}(q^{2})$ is
more flat in the whole $q^{2}$ region.

  * •

Our predictions of $f^{\eta}_{+}(q^{2})$ are in good agreement with the recent
BESIII predictions for the $\eta\to\gamma\gamma$ channel.

  * •

In the LCSR $q^{2}$-region, our results have good agreement with the LCSR 2013
and 2015, the CCQM, and the two sets of BES-III predictions within errors.

  * •

The SSE of $f^{\eta}_{+}(q^{2})$ in the region of $q^{2}\in[1.2,2.0]~{}{\rm
GeV}^{2}$ have agreement with the LCSR in 2015 results within errors. However,
our predictions are larger than the LCSR in 2013 and CCQM predictions due to
the different $\eta$-meson distribution amplitudes or different method.

Figure 8: Decay width for the
$D_{s}^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}(\ell=e,\mu)$ versus $q^{2}$
within uncertainties. The LCSR in 2013 Offen:2013nma  and 2015
Duplancic:2015zna , the CCQM Ivanov:2019nqd , and two sets of BESIII
collaboration Ablikim:2019rjz  predictions are also present as a comparison.

###  III.4 Decay widthes and branching fractions for the semileptonic decay
$D_{s}^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$

One can get the differential decay widths by using the formula (56). For the
CKM matrix element $|V_{cs}|$, we take it to be the average value of leptonic
and semileptonic decay $c\to s$ processes coming from PDG Narison:2014ska ,
i.e. $|V_{cs}|=0.987\pm 0.011$. After taking the derived
$D_{s}\to{\eta^{(\prime)}}$ TFFs into the decay widths, we present the
differential decay widths in Figure 8. As a comparison, we also give the LCSR
in 2013 Offen:2013nma  and 2015 Duplancic:2015zna , the CCQM Ivanov:2019nqd ,
and two sets of BESIII collaboration Ablikim:2019rjz  predictions. The LCSR
and the CCQM results are calculated by applying their TFFs into the width
formula. Those figures show that our prediction for
$D_{s}^{+}\to\eta\ell^{+}\nu_{\ell}$ is in agreement with LCSR 2015, CCQM, the
BESIII-I and the BESIII-II results within errors, and the
$D_{s}^{+}\to\eta^{\prime}\ell^{+}\nu_{\ell}$ agrees with the LCSR 2013 and
2015, the CCQM, the BESIII-I and BESIII-II predictions within errors. All the
results are convergence to zero at the small recoil region
$q^{2}_{\max}=(m_{D_{s}}-m_{\eta^{(\prime)}})^{2}$, which indicates that our
results are reasonable.

Table 6: Branching factions of $D_{s}^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$ with $\ell=e$ and $\mu$ (in unit $10^{-2}$). The errors are squared averages of all the mentioned error sources. As a comparison, we also present the predictions for various methods. | Mode | ${\cal B}(D_{s}^{+}\to\eta e^{+}\nu_{e})$ | ${\cal B}(D_{s}^{+}\to\eta\mu^{+}\nu_{\mu})$
---|---|---|---
| BESIII Ablikim:2019rjz ; Ablikim:2017omq  | $2.323\pm 0.063\pm 0.063$ | $2.42\pm 0.46\pm 0.11$
| CLEO Hietala:2015jqa  | $2.28\pm 0.14\pm 0.19$ | -
Experimental results | CLEO Yelton:2009aa  | $2.48\pm 0.29\pm 0.13$ | -
| PDG Zyla:2020zbs  | $2.32\pm 0.08$ | $2.4\pm 0.5$
| This work (LCSR) | $2.346_{-0.331}^{+0.418}$ | $2.320_{-0.327}^{+0.413}$
| LFQM Cheng:2017pcq  | $2.26\pm 0.21$ | $2.22\pm 0.20$
| CCQM Ivanov:2019nqd  | $2.24$ | $2.18$
Theoretical predictions | LCSR Offen:2013nma  | $2.00\pm 0.32$ | -
| LCSR Duplancic:2015zna  | $2.40\pm 0.28$ | -
| QCD SR-I Colangelo:2001cv  | $2.6\pm 0.7$ | -
| QCD SR-II Colangelo:2001cv  | $2.3\pm 0.4$ | -
| Mode | ${\cal B}(D_{s}^{+}\to\eta^{\prime}e^{+}\nu_{e})$ | ${\cal B}(D_{s}^{+}\to\eta^{\prime}\mu^{+}\nu_{\mu})$
| BESIII Ablikim:2019rjz ; Ablikim:2017omq  | $0.824\pm 0.073\pm 0.027$ | $1.06\pm 0.54\pm 0.07$
| CLEO Hietala:2015jqa  | $0.68\pm 0.15\pm 0.06$ | -
Experimental results | CLEO Yelton:2009aa  | $0.91\pm 0.33\pm 0.05$ | -
| PDG Zyla:2020zbs  | $0.80\pm 0.07$ | $1.1\pm 0.5$
| This work (LCSR) | $0.792_{-0.118}^{+0.141}$ | $0.773_{-0.115}^{+0.138}$
| LFQM Cheng:2017pcq  | $0.89\pm 0.09$ | $0.85\pm 0.08$
| CCQM Ivanov:2019nqd  | $0.83$ | $0.79$
Theoretical predictions | LCSR Offen:2013nma  | $0.75\pm 0.23$ | -
| LCSR Duplancic:2015zna  | $0.79\pm 0.14$ | -
| QCD SR-I Colangelo:2001cv  | $0.89\pm 0.34$ | -
| QCD SR-II Colangelo:2001cv  | $1.0\pm 0.2$ | -

After integrating over the whole $m_{\ell}^{2}\leq
q^{2}\leq(m_{D_{s}}-m_{\eta^{(\prime)}})^{2}$ region for the differential
decay widths, we obtain the total decay widths for
$D_{s}^{+}\to\eta^{(\prime)}\ell^{+}\nu_{\ell}$ with two different channel
$\Gamma(D_{s}^{+}\to{\eta^{(\prime)}}e^{+}\nu_{e})$ and
$\Gamma(D_{s}^{+}\to{\eta^{(\prime)}}\mu^{+}\nu_{\mu})$, i.e.,

$\displaystyle\Gamma(D_{s}^{+}\to\eta
e^{+}\nu_{e})=30.634_{-4.323}^{+5.453}\times 10^{-15}~{}{\rm GeV},$
$\displaystyle\Gamma(D_{s}^{+}\to\eta\mu^{+}\nu_{\mu})=30.298_{-4.275}^{+5.395}\times
10^{-15}~{}{\rm GeV},$
$\displaystyle\Gamma(D_{s}^{+}\to\eta^{\prime}e^{+}\nu_{e})=10.345_{-1.534}^{+1.844}\times
10^{-15}~{}{\rm GeV},$
$\displaystyle\Gamma(D_{s}^{+}\to\eta^{\prime}\mu^{+}\nu_{\mu})=10.096_{-1.498}^{+1.800}\times
10^{-15}~{}{\rm GeV}.$ (76)

Then, by using the lifetime of the initial state $D_{s}^{+}$-meson,
$\tau_{D_{s}^{+}}=(0.504\pm 0.007)~{}{\rm ps}$ Zyla:2020zbs , the branching
fractions for the two different semileptonic decay channels
$D_{s}\to{\eta^{(\prime)}}{{\ell}^{+}}\nu_{\ell}$ with $\ell=(e,\mu)$ can be
obtained, which are presented in Table 6. Here, we also listed the BESIII
Ablikim:2019rjz ; Ablikim:2017omq , the PDG Zyla:2020zbs , the CLEO
Yelton:2009aa  for the experimental results, and the CCQM Ivanov:2019nqd , the
LFQM Cheng:2017pcq , the QCDSR-I, II Colangelo:2001cv , the LCSR
Duplancic:2015zna ; Offen:2013nma  for theoretical predictions. Our results
are closer to the BESIII, PDG, CLEO results, all of which are within $1\sigma$
uncertainties.

Figure 9: $f^{\eta}_{+}(0)/f^{\eta^{\prime}}_{+}(0)$ as a function of the Borel parameter $M^{2}$, where the shaded band is induced by the variations of squared average of all input parameters Table 7: The $\varphi$ with respect to mixing angle, ${\cal R}_{\eta^{\prime}/\eta}^{\ell}$ for different models and experimental values. As a comparison, we also present the experimental and theoretical predictions. Mode | Angle $\varphi$ | ${\cal R}_{\eta^{\prime}/\eta}^{\ell}$
---|---|---
This work (LCSR) with $\ell=e$ | ${41.2^{\circ}}_{-0.06}^{+0.05}$ | $0.338_{-0.051}^{+0.057}$
This work (LCSR) with $\ell=\mu$ | ${41.2^{\circ}}_{-0.06}^{+0.05}$ | $0.333_{-0.058}^{+0.050}$
CLEO Brandenburg:1995qq  | - | $0.35\pm 0.09\pm 0.07$
BESIII Ablikim:2016rqq  | - | $0.40\pm 0.14\pm 0.02$
LFQM Wei:2009nc  | ${39}^{\circ}$ | $0.39$
LCSR Duplancic:2015zna  | ${41.8}^{\circ}$ | $0.33\pm 0.07$
QCD SR Colangelo:2001cv  | ${40}^{\circ}$ | $0.44\pm 0.01$
KLOE Ambrosino:2006gk  | ${41.4}^{\circ}$ | -
LCSR-I Azizi:2010zj  | ${39.7}^{\circ}$ | $0.32\pm 0.02$
LCSR-II Azizi:2010zj  | ${41.5}^{\circ}$ | $0.27\pm 0.01$

After substituting the corresponding terms of Eq. (9), one can get the mixing
angle of $\eta-\eta^{\prime}$. The mixing angle $\tan\varphi$ with the
parameter of Borel parameter $M^{2}$ are shown in Figure 9. In the whole Borel
parameter $M^{2}$ region, the mixing angle is changed slightly, which also
indicates that there have a stable Borel window for $\tan\varphi$. Numerical
results for the single mixing angle $\varphi$ are present in Table 7. Our
result is closer to the KLOE Ambrosino:2006gk  and the LCSR predictions
Duplancic:2015zna ; Azizi:2010zj , which is slightly larger than other LFQM
Wei:2009nc  and QCD SR Colangelo:2001cv  predictions.

Furthermore, it is useful to study the ratio for the different decay channel
${\cal R}_{\eta^{\prime}/\eta}^{\ell}$ related to the mixing angle, which has
the basic definition

$\displaystyle{\cal R}_{\eta^{\prime}/\eta}^{\ell}=\frac{{\cal
B}(D_{s}\to\eta^{\prime}\ell^{+}\nu_{\ell})}{{\cal
B}(D_{s}\to\eta\ell^{+}\nu_{\ell})}.$ (77)

Numerical results together with different experimental and theoretical
predictions are given in Table 7. Our results are in agreement with the CLEO
Brandenburg:1995qq  and the BESIII Ablikim:2016rqq  predictions, and the LCSR
predictions Duplancic:2015zna ; Azizi:2010zj . This can be considered as a
good test of the correctness of the considered internal structure for the
$D_{s}$-meson as well as the mixing angle between $\eta$ and $\eta^{\prime}$
states.

##  IV Summary

In this paper, we have calculated the moments
$\langle\xi_{2;\eta^{(\prime)}}^{n}\rangle|_{\mu}$ of $\eta^{(\prime)}$-meson
LCDA with $n=(2,4,6)$ up to NLO correction and completely dimension-six
condensates within BFTSR, which are shown in Eq. (39). Due to the
$\eta^{(\prime)}$-meson should be considered as $s\bar{s}$ component in
$D_{s}\to\eta^{(\prime)}$ decay processes, we have also taken the
$I_{m_{s}^{2}}(n,M^{2})$ corrections into consideration, i.e. Eq. (40), the
detailed terms of OPE are listed in the Appendix B.

Then, we have sought a reasonable continuum threshold $s_{0}$ and stable Borel
windows for $\langle\xi_{2;\eta^{(\prime)}}^{n}\rangle|_{\mu}$ with
$n=(2,4,6)$ by using the traditional three criteria for the SVZ sum rules,
which are present in Figs. 1, 2 and Table 1. By using the expression between
two different Gegenbauer and LCDA moments, i.e. Eq. (43), we have presented
the first three Gegenbaner and LCDA moments
$a^{n}_{2;\eta^{(\prime)}}(\mu_{0})$ and
$\langle\xi^{n}_{2;\eta^{(\prime)}}\rangle|_{\mu_{0}}$ with $n=(2,4,6)$ for
the $\phi_{2;\eta^{(\prime)}}(u,\mu_{0})$ within errors in Table 2. Our
results are in agreement with the CLEO fit and the BABAR fit predictions.
Meanwhile, we have exhibited the curves of
$\phi_{2;\eta^{(\prime)}}(u,\mu_{0})$ of our prediction with $n=(2,4,6)$ and
compare with others.

Furthermore, the TFFs $f^{\eta^{(\prime)}}_{+}(q^{2})$ have been given in Eq.
(LABEL:Eq:fp) up to NLO QCD corrections for twist-2, 3 LCDA contributions. The
TFFs at large recoil region have been presented in Table 3 with respect to
other theoretical predictions. After extrapolating it to the whole physical
$q^{2}$-region via simplified series expansion, we have shown the behavior of
TFFs in Figure 7. The differential/total decay widths and branching fractions
in this work have also been given. Our results are in agreement with the
BESIII and the PDG average value within errors. Finally, we have presented the
mixing angle $\tan\varphi$ and ratio for different decay channel ${\cal
R}_{\eta^{\prime}/\eta}^{\ell}$, which agree with theoretical and experimental
results within errors. Thus, the QCDSR within BFTSR can be considered a good
approach in dealing with the heavy-to-light semileptonic processes, and we
hope more data can be achieved in the near future for more precise studies.

##  V Acknowledgments

We are grateful to Dr. Xu-Chang Zheng for helpful discussions and valuable
suggestions. Hai-Bing Fu would like to thank the Institute of Theoretical
Physics in Chongqing University (CQUITP) for kind hospitality. This work was
supported in part by the National Natural Science Foundation of China under
Grant No.11765007, No.11947406, No.11625520, and No.12047564, the Project of
Guizhou Provincial Department of Science and Technology under Grant
No.KY[2019]1171, and No.ZK[2021]024, the Project of Guizhou Provincial
Department of Education under Grant No.KY[2021]030 and No.KY[2021]003, the
China Postdoctoral Science Foundation under Grant Nos.2019TQ0329, 2020M670476,
the Chongqing Graduate Research and Innovation Foundation under Grant
No.ydstd1912, the Fundamental Research Funds for the Central Universities
under Grant No.2020CQJQY-Z003, and the Project of Guizhou Minzu University
under Grant No. GZMU[2019]YB19.

##  Appendix A Quark propagators and vertex for background field theory
framework

###  A.1 Quark propagators

Quarks field and its propagators satisfy the equation Eq. (17) in the BFT,
slove to Eq. (17), we have the following form

$\displaystyle S_{F}(x,0)=(i/\\!\\!\\!\\!D+m){\rm\cal D}(x,0),$ (78)

where ${\rm\cal D}(x,0)$ satisfies the equation

$\displaystyle(\partial^{2}-{\cal P}_{\mu}\partial^{\mu}-{\cal Q}+m^{2}){\cal
D}(x,0)=\delta^{4}(x),$ (79)

with ${\cal P}_{\mu}=2i{\cal A}_{\mu}$, ${\cal
Q}=\gamma^{\nu}\gamma^{\mu}[{\cal A}_{\nu}(x){\cal
A}_{\mu}(x)+i\partial_{\nu}{\cal A}_{\mu}(x)]$. After using above equation up
here we can get quark propagator $S_{F}(x,0)$ subexpression

$\displaystyle
S_{F}(x,0)=S_{F}^{0}(x,0)+S_{F}^{2}(x,0)+S_{F}^{3}(x,0)+\sum\limits_{i=1}^{2}S_{F}^{4(i)}(x,0)+\sum\limits_{i=1}^{3}S_{F}^{5(i)}(x,0)+\sum\limits_{i=1}^{5}S_{F}^{6(i)}(x,0).$
(80)

In order to have a clear look at the detailed quark propagator $S_{F}(x,0)$
separately, we list the expression of $S_{F}^{d}(x,0)$ with $d\leq 6$ as
follows,

Figure 10: Feynman diagrams for the quark propagator within the background
field theory which shall results in operators up to dimension-six. The symbol
“$\times$” attached to the gluon line indicates the tensor of the local gluon
background field, in which “$n$” stands for $n$-th order covariant derivative.

$\displaystyle S_{F}^{0}(x,0)=i\int\frac{d^{4}p}{(2\pi)^{4}}e^{-ip\cdot
x}\bigg{[}-\frac{m+/\\!\\!\\!p}{m^{2}-p^{2}}\bigg{]},$ (81) $\displaystyle
S_{F}^{2}(x,0)=i\int\frac{d^{4}p}{(2\pi)^{4}}e^{-ip\cdot
x}\bigg{[}-\frac{i}{2}\frac{\gamma^{\mu}(m-/\\!\\!\\!p)\gamma^{\nu}}{(m^{2}-p^{2})^{2}}G_{\mu\nu}\bigg{]},$
(82) $\displaystyle S_{F}^{3}(x,0)=i\int\frac{d^{4}p}{(2\pi)^{4}}e^{-ip\cdot
x}\bigg{\\{}-\frac{2}{3}\bigg{[}\frac{(\gamma^{\mu}p^{\rho}+\gamma^{\rho}p^{\mu})(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{3}}-\frac{g^{\mu\rho}}{(m^{2}-p^{2})^{2}}\bigg{]}\gamma^{\nu}G_{\mu\nu;\rho}\bigg{\\}},$
(83) $\displaystyle
S_{F}^{4(1)}(x,0)=i\int\frac{d^{4}p}{(2\pi)^{4}}e^{-ip\cdot
x}\bigg{\\{}\frac{1}{4}\bigg{[}\frac{\gamma^{\mu}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{3}}-\frac{2p^{\mu}}{(m^{2}-p^{2})^{3}}\bigg{]}\gamma^{\nu}\gamma^{\rho}\gamma^{\sigma}+\frac{1}{2}~{}\bigg{[}\frac{(m+/\\!\\!\\!p)\gamma^{\mu}}{(m^{2}-p^{2})^{3}}~{}g^{\nu\sigma}$
$\displaystyle\hskip
48.36958pt+4\frac{\gamma^{\mu}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{4}}p^{\nu}p^{\sigma}\bigg{]}\gamma^{\rho}\bigg{\\}}G_{\mu\nu}G_{\rho\sigma},$
(84) $\displaystyle
S_{F}^{4(2)}(x,0)=i\int\frac{d^{4}p}{(2\pi)^{4}}e^{-ip\cdot
x}\bigg{\\{}~{}\frac{i}{4}\bigg{[}\frac{{g^{\\{\mu\rho}}\gamma^{\sigma\\}}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{3}}-\frac{2g^{\\{\mu\rho}p^{\sigma\\}}}{(m^{2}-p^{2})^{3}}+4\frac{{\gamma^{\\{\mu}}p^{\rho}p^{\sigma\\}}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{4}}\bigg{]}\gamma^{\nu}$
$\displaystyle\hskip 48.36958pt\times G_{\mu\nu;\rho\sigma}\bigg{\\}},$ (85)
$\displaystyle S_{F}^{5(1)}(x,0)=i\int\frac{d^{4}p}{(2\pi)^{4}}e^{-ip\cdot
x}\bigg{\\{}\\!-\frac{i}{3}\bigg{[}\bigg{(}3\frac{\gamma^{\mu}(m-/\\!\\!\\!p)\gamma^{\nu}}{(m^{2}-p^{2})^{4}}(p^{\lambda}\gamma^{\rho}+p^{\rho}\gamma^{\lambda})-\frac{\gamma^{\nu}({g^{\mu\lambda}}\gamma^{\rho}+g^{\mu\rho}\gamma^{\lambda})}{(m^{2}-p^{2})^{3}}\bigg{)}$
$\displaystyle\hskip
48.36958pt\cdot~{}\gamma^{\sigma}+4~{}\bigg{(}\frac{{\gamma^{\mu}}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{4}}{g^{\\{\nu\sigma}}{p^{\lambda\\}}}+2\frac{p^{\mu}{g^{\\{\nu\sigma}}p^{\lambda\\}}}{(m^{2}-p^{2})^{4}}+6~{}\frac{\gamma^{\mu}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{5}}p^{\nu}p^{\sigma}p^{\lambda}\bigg{)}\gamma^{\rho}\bigg{]}~{}G_{\mu\nu}$
$\displaystyle\hskip 48.36958pt\times G_{\rho\sigma;\lambda}\bigg{\\}},$ (86)
$\displaystyle S_{F}^{5(2)}(x,0)=i\int{\frac{d^{4}p}{(2\pi)^{4}}}e^{-ip\cdot
x}\bigg{\\{}\frac{2i}{3}\bigg{[}\bigg{(}\frac{g^{\mu\lambda}}{(m^{2}-p^{2})^{3}}+6\frac{p^{\mu}p^{\lambda}}{(m^{2}-p^{2})^{4}}\bigg{)}\gamma^{\nu}\gamma^{\rho}\gamma^{\sigma}-2~{}\bigg{(}\frac{\gamma^{\mu}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{4}}$
$\displaystyle\hskip
48.36958pt\times{g^{\\{\nu\sigma}}{p^{\lambda\\}}}+2\frac{p^{\mu}{g^{\\{\nu\sigma}}p^{\lambda\\}}}{(m^{2}-p^{2})^{4}}+6\frac{\gamma^{\mu}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{5}}p^{\nu}p^{\sigma}p^{\lambda}\bigg{)}\gamma^{\rho}\bigg{]}G_{\mu\nu;\lambda}G_{\rho\sigma}\bigg{\\}},$
(87) $\displaystyle
S_{F}^{5(3)}(x,0)=i\int{\frac{d^{4}p}{(2\pi)^{4}}}e^{-ip\cdot
x}~{}\frac{4}{15}\bigg{[}\frac{{g^{\\{\rho\sigma}}p^{\lambda}\gamma^{\mu\\}}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{4}}-\frac{2{g^{\\{\rho\sigma}}p^{\lambda}p^{\mu\\}}}{(m^{2}-p^{2})^{4}}+6\frac{{\gamma^{\\{\mu}}p^{\rho}p^{\sigma}p^{\lambda\\}}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{5}}$
$\displaystyle\hskip
48.36958pt-\frac{g^{(\mu\nu\sigma\lambda)}}{(m^{2}-p^{2})^{3}}\bigg{]}\gamma^{\nu}G_{\mu\nu;\rho\sigma\lambda},$
(88) $\displaystyle
S_{F}^{6(1)}(x,0)=i\int{\frac{d^{4}p}{(2\pi)^{4}}}e^{-ip\cdot
x}\frac{i}{8}\bigg{\\{}\bigg{[}\frac{\gamma^{\mu}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{4}}-\\!4\frac{p^{\mu}}{(m^{2}-p^{2})^{4}}\bigg{]}\gamma^{\nu}\gamma^{\rho}\gamma^{\sigma}\gamma^{\lambda}\gamma^{\tau}\\!+\\!2\bigg{[}\frac{3\gamma^{\mu}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{4}}$
$\displaystyle\hskip 48.36958pt\times
g^{\sigma\tau}+16\frac{\gamma^{\mu}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{5}}p^{\sigma}p^{\tau}\\!-\\!4\frac{g^{\mu\sigma}p^{\tau}+g^{\mu\tau}p^{\sigma}}{(m^{2}-p^{2})^{4}}\bigg{]}\gamma^{\nu}\gamma^{\rho}\gamma^{\lambda}\bigg{\\}}G_{\mu\nu}G_{\rho\sigma}G_{\lambda\tau},$
(89) $\displaystyle
S_{F}^{6(2)}(x,0)=i\int{\frac{d^{4}p}{(2\pi)^{4}}}e^{-ip\cdot
x}\bigg{(}-\frac{1}{8}\bigg{)}\bigg{\\{}\bigg{[}3\frac{\gamma^{\mu}(m-/\\!\\!\\!p)\gamma^{\nu}}{(m^{2}-p^{2})^{4}}{g^{\\{\lambda\tau}}\gamma^{\rho\\}}+16\frac{\gamma^{\mu}(m-/\\!\\!\\!p)\gamma^{\nu}}{(m^{2}-p^{2})^{5}}~{}\gamma^{\\{\rho}p^{\lambda}p^{\tau\\}}$
$\displaystyle\hskip
48.36958pt-4\frac{\gamma^{\nu}}{(m^{2}-p^{2})^{4}}{g^{\mu\\{\lambda}}p^{\tau}\gamma^{\rho\\}}\bigg{]}{\gamma^{\sigma}}+4\bigg{[}\frac{m+/\\!\\!\\!p}{(m^{2}-p^{2})^{4}}~{}g^{(\nu\sigma\tau\lambda)}+6\frac{m+/\\!\\!\\!p}{(m^{2}-p^{2})^{5}}~{}\,{g^{\\{\nu\sigma}}p^{\tau}p^{\lambda\\}}$
$\displaystyle\hskip
48.36958pt+48\frac{m+/\\!\\!\\!p}{(m^{2}-p^{2})^{6}}p^{\nu}p^{\sigma}p^{\tau}p^{\lambda}\bigg{]}\gamma^{\mu}\gamma^{\rho}\bigg{\\}}G_{\mu\nu}G_{\rho\sigma;\lambda\tau},$
(90) $\displaystyle
S_{F}^{6(3)}(x,0)=i\int{\frac{d^{4}p}{(2\pi)^{4}}}e^{-ip\cdot
x}\bigg{(}-\frac{2}{9}\bigg{)}\bigg{\\{}3\bigg{[}\bigg{(}2\frac{{\gamma^{\mu}}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{5}}p^{\lambda}p^{\tau}-\frac{{g^{\\{\mu\lambda}}p^{\tau\\}}}{(m^{2}-p^{2})^{4}}-\frac{4p^{\mu}p^{\lambda}p^{\tau}}{(m^{2}-p^{2})^{5}}\bigg{)}$
$\displaystyle\hskip
48.36958pt\cdot~{}\gamma^{\nu}\gamma^{\rho}+(\mu\leftrightarrow\lambda)+(\rho\leftrightarrow\tau)+(\mu\leftrightarrow\lambda,\rho\leftrightarrow\tau)\bigg{]}\gamma^{\sigma}+4\,\bigg{[}\frac{m+/\\!\\!\\!p}{(m^{2}-p^{2})^{4}}g^{(\nu\lambda\sigma\tau)}+6$
$\displaystyle\hskip
48.36958pt\times\frac{m+/\\!\\!\\!p}{(m^{2}-p^{2})^{5}}g^{\\{\nu\lambda}{p^{\sigma}}p^{\tau\\}}+48\frac{m+/\\!\\!\\!p}{(m^{2}-p^{2})^{6}}p^{\nu}p^{\lambda}p^{\sigma}p^{\tau}\bigg{]}\gamma^{\mu}\gamma^{\rho}\bigg{\\}}G_{\mu\nu;\lambda}G_{\rho\sigma;\tau},$
(91) $\displaystyle
S_{F}^{6(4)}(x,0)=i\int{\frac{d^{4}p}{(2\pi)^{4}}}e^{-ip\cdot
x}~{}\bigg{(}-\frac{1}{2}\bigg{)}\bigg{\\{}\bigg{[}\frac{m+/\\!\\!\\!p}{(m^{2}-p^{2})^{4}}g^{\\{\nu\tau\lambda\sigma\\}}+6\frac{m+/\\!\\!\\!p}{(m^{2}-p^{2})^{5}}{g^{\\{\nu\tau}}p^{\lambda}p^{\sigma\\}}+48$
$\displaystyle\hskip
48.36958pt\times\frac{(m+/\\!\\!\\!p)p^{\nu}p^{\tau}p^{\lambda}p^{\sigma}}{(m^{2}-p^{2})^{6}}\bigg{]}\gamma^{\mu}\gamma^{\rho}-3\bigg{[}\frac{{g^{\\{\mu\lambda}}p^{\tau\\}}}{(m^{2}-p^{2})^{4}}+\frac{8p^{\mu}p^{\lambda}p^{\tau}}{(m^{2}-p^{2})^{5}}\bigg{]}\gamma^{\nu}\gamma^{\rho}\gamma^{\sigma}\bigg{\\}}G_{\mu\nu;\lambda\tau}G_{\rho\sigma},$
(92) $\displaystyle
S_{F}^{6(5)}(x,0)=i\int{\frac{d^{4}p}{(2\pi)^{4}}}e^{-ip\cdot
x}\bigg{\\{}-\frac{i}{18}~{}\bigg{[}\frac{{g^{[\rho\sigma\lambda\tau}}{\gamma^{\mu]}}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{4}}~{}-4~{}\frac{{g^{\\{\rho\sigma\lambda\tau}}p^{\mu\\}}}{(m^{2}-p^{2})^{4}}+6~{}g^{\\{\rho\sigma}p^{\lambda}p^{\tau}\gamma^{\mu\\}}$
$\displaystyle\hskip
48.36958pt\times\frac{m-/\\!\\!\\!p}{(m^{2}-p^{2})^{5}}-12\frac{{g^{\\{\rho\sigma}}p^{\lambda}p^{\tau}p^{\mu\\}}}{(m^{2}-p^{2})^{5}}+48\frac{{p^{\\{\rho}}p^{\sigma}p^{\lambda}p^{\tau}p^{\mu\\}}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{6}}\bigg{]}\gamma^{\nu}G_{\mu\nu;\rho\sigma\lambda\tau}\bigg{\\}}.$
(93)

The Feynman diagrams for the quark propagators Eqs. (81)-(93) that with
various gauge invariant tensors are shown in Figure 10, where thirteen
figures, i.e. Figure 10a, $...$, 10m, correspond to $S_{F}^{0}(x,0)$, $...$,
$S_{F}^{6(5)}(x,0)$, respectively. The symbol “$\times$” attached to the gluon
line indicates the tensor of the local gluon background field with “$n$”
stands for $n$-th order covariant derivative. For quark propagator from $x$ to
$0$, the following relation can be used,

$\displaystyle S_{F}(0,x|{\cal A})=CS_{F}^{T}(x,0|-{\cal A}^{T}){C^{-1}}.$
(94)

The detailed expression are

$\displaystyle S_{F}^{0}(0,x)=i\int\frac{d^{4}p}{(2\pi)^{4}}e^{-ip\cdot
x}\bigg{[}-\frac{m-/\\!\\!\\!p}{m^{2}-p^{2}}\bigg{]},$ (95) $\displaystyle
S_{F}^{2}(0,x)=i\int\frac{d^{4}p}{(2\pi)^{4}}e^{-ip\cdot
x}\bigg{[}\frac{i}{2}\frac{\gamma^{\nu}(m+/\\!\\!\\!p)\gamma^{\mu}}{(m^{2}-p^{2})^{2}}G_{\mu\nu}\bigg{]},$
(96) $\displaystyle S_{F}^{3}(0,x)=i\int\frac{d^{4}p}{(2\pi)^{4}}e^{-ip\cdot
x}\bigg{\\{}\frac{2}{3}\gamma^{\nu}\bigg{[}\frac{(m+/\\!\\!\\!p)(\gamma^{\mu}p^{\rho}+\gamma^{\rho}p^{\mu})}{(m^{2}-p^{2})^{3}}+\frac{g^{\mu\rho}}{(m^{2}-p^{2})^{2}}\bigg{]}G_{\mu\nu;\rho}\bigg{\\}},$
(97) $\displaystyle
S_{F}^{4(1)}(0,x)=i\int\frac{d^{4}p}{(2\pi)^{4}}e^{-ip\cdot
x}\bigg{\\{}\frac{1}{4}\gamma^{\sigma}\gamma^{\rho}\gamma^{\nu}\bigg{[}\frac{(m+/\\!\\!\\!p)\gamma^{\mu}}{(m^{2}-p^{2})^{3}}+\frac{2p^{\mu}}{(m^{2}-p^{2})^{3}}\bigg{]}+\frac{1}{2}\gamma^{\rho}\bigg{[}\frac{\gamma^{\mu}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{3}}~{}g^{\nu\sigma}$
$\displaystyle\hskip
48.36958pt+4\frac{(m+/\\!\\!\\!p)\gamma^{\mu}}{(m^{2}-p^{2})^{4}}p^{\nu}p^{\sigma}\bigg{]}\bigg{\\}}G_{\rho\sigma}G_{\mu\nu},$
(98) $\displaystyle
S_{F}^{4(2)}(0,x)=i\int\frac{d^{4}p}{(2\pi)^{4}}e^{-ip\cdot
x}\bigg{\\{}-\frac{i}{4}\gamma^{\nu}\bigg{[}\frac{(m+/\\!\\!\\!p){g^{\\{\mu\rho}}\gamma^{\sigma\\}}}{(m^{2}-p^{2})^{3}}+\frac{2g^{\\{\mu\rho}p^{\sigma\\}}}{(m^{2}-p^{2})^{3}}+4\frac{(m+/\\!\\!\\!p){\gamma^{\\{\mu}}p^{\rho}p^{\sigma\\}}}{(m^{2}-p^{2})^{4}}\bigg{]}$
$\displaystyle\hskip 48.36958pt\times G_{\mu\nu;\rho\sigma}\bigg{\\}},$ (99)
$\displaystyle S_{F}^{5(1)}(0,x)=i\int\frac{d^{4}p}{(2\pi)^{4}}e^{-ip\cdot
x}\bigg{\\{}\\!-\frac{i}{3}\bigg{[}\gamma^{\sigma}\bigg{(}-3(p^{\lambda}\gamma^{\rho}+p^{\rho}\gamma^{\lambda})\frac{\gamma^{\nu}(m+/\\!\\!\\!p)\gamma^{\mu}}{(m^{2}-p^{2})^{4}}-\frac{({g^{\mu\lambda}}\gamma^{\rho}+g^{\mu\rho}\gamma^{\lambda})}{(m^{2}-p^{2})^{3}}$
$\displaystyle\hskip
48.36958pt~{}\cdot~{}\gamma^{\nu}\bigg{)}+4\gamma^{\rho}\bigg{(}\frac{{\gamma^{\mu}}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{4}}{g^{\\{\nu\sigma}}{p^{\lambda\\}}}+2\frac{p^{\mu}{g^{\\{\nu\sigma}}p^{\lambda\\}}}{(m^{2}-p^{2})^{4}}+6\frac{\gamma^{\mu}(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{5}}p^{\nu}p^{\sigma}p^{\lambda}\bigg{)}\bigg{]}G_{\rho\sigma;\lambda}$
$\displaystyle\hskip 48.36958pt\times G_{\mu\nu}\bigg{\\}},$ (100)
$\displaystyle S_{F}^{5(2)}(0,x)=i\int{\frac{d^{4}p}{(2\pi)^{4}}}e^{-ip\cdot
x}\frac{2i}{3}\bigg{[}\gamma^{\sigma}\gamma^{\rho}\gamma^{\nu}\bigg{(}\frac{g^{\mu\lambda}}{(m^{2}-p^{2})^{3}}+6\frac{p^{\mu}p^{\lambda}}{(m^{2}-p^{2})^{4}}\bigg{)}\\!-2\gamma^{\rho}\bigg{(}\\!\\!-\frac{(m+/\\!\\!\\!p)\gamma^{\mu}}{(m^{2}-p^{2})^{4}}$
$\displaystyle\hskip
48.36958pt\times{g^{\\{\nu\sigma}}{p^{\lambda\\}}}-2\frac{p^{\mu}{g^{\\{\nu\sigma}}p^{\lambda\\}}}{(m^{2}-p^{2})^{4}}-6\frac{(m+/\\!\\!\\!p)\gamma^{\mu}}{(m^{2}-p^{2})^{5}}p^{\nu}p^{\sigma}p^{\lambda}\bigg{)}\bigg{]}G_{\rho\sigma}G_{\mu\nu;\lambda},$
(101) $\displaystyle
S_{F}^{5(3)}(0,x)=i\int{\frac{d^{4}p}{(2\pi)^{4}}}e^{-ip\cdot
x}~{}\frac{4}{15}~{}\gamma^{\nu}~{}\bigg{[}-\frac{(m+/\\!\\!\\!p){g^{\\{\rho\sigma}}p^{\lambda}\gamma^{\mu\\}}}{(m^{2}-p^{2})^{4}}~{}-~{}\frac{2{g^{\\{\rho\sigma}}p^{\lambda}p^{\mu\\}}}{(m^{2}-p^{2})^{4}}-6~{}\frac{(m-/\\!\\!\\!p)}{(m^{2}-p^{2})^{5}}$
$\displaystyle\hskip
48.36958pt~{}\cdot~{}{\gamma^{\\{\mu}}p^{\rho}p^{\sigma}p^{\lambda\\}}-\frac{g^{(\mu\nu\sigma\lambda)}}{(m^{2}-p^{2})^{3}}\bigg{]}G_{\mu\nu;\rho\sigma\lambda},$
(102) $\displaystyle
S_{F}^{6(1)}(0,x)=i\int{\frac{d^{4}p}{(2\pi)^{4}}}e^{-ip\cdot
x}\Big{(}-\frac{i}{8}\Big{)}\bigg{\\{}\gamma^{\tau}\gamma^{\lambda}\gamma^{\sigma}\gamma^{\rho}\gamma^{\nu}\bigg{[}\frac{(m+/\\!\\!\\!p)\gamma^{\mu}}{(m^{2}-p^{2})^{4}}+4\frac{p^{\mu}}{(m^{2}-p^{2})^{4}}\bigg{]}+2\gamma^{\lambda}\gamma^{\rho}\gamma^{\nu}$
$\displaystyle\hskip
48.36958pt\times\bigg{[}\frac{3(m+/\\!\\!\\!p)\gamma^{\mu}}{(m^{2}-p^{2})^{4}}g^{\sigma\tau}+16\frac{(m+/\\!\\!\\!p)\gamma^{\mu}}{(m^{2}-p^{2})^{5}}p^{\sigma}p^{\tau}\\!+\\!4\frac{g^{\mu\sigma}p^{\tau}+g^{\mu\tau}p^{\sigma}}{(m^{2}-p^{2})^{4}}\bigg{]}\bigg{\\}}G_{\lambda\tau}G_{\rho\sigma}G_{\mu\nu},$
(103) $\displaystyle
S_{F}^{6(2)}(0,x)=i\int{\frac{d^{4}p}{(2\pi)^{4}}}e^{-ip\cdot
x}\bigg{(}-\frac{1}{8}\bigg{)}\bigg{\\{}\gamma^{\sigma}\\!\bigg{[}3{g^{\\{\lambda\tau}}\gamma^{\rho\\}}\frac{\gamma^{\nu}(m+/\\!\\!\\!p)\gamma^{\mu}}{(m^{2}-p^{2})^{4}}+16\gamma^{\\{\rho}p^{\lambda}p^{\tau\\}}\frac{\gamma^{\nu}(m+/\\!\\!\\!p)\gamma^{\mu}}{(m^{2}-p^{2})^{5}}$
$\displaystyle\hskip
48.36958pt+4\frac{{g^{\mu\\{\lambda}}p^{\tau}\gamma^{\rho\\}}\gamma^{\nu}}{(m^{2}-p^{2})^{4}}\bigg{]}+4~{}\gamma^{\rho}\gamma^{\mu}~{}\bigg{[}\frac{m-/\\!\\!\\!p}{(m^{2}-p^{2})^{4}}~{}g^{(\nu\sigma\tau\lambda)}~{}+6\frac{m-/\\!\\!\\!p}{(m^{2}-p^{2})^{5}}~{}{g^{\\{\nu\sigma}}p^{\tau}p^{\lambda\\}}+48$
$\displaystyle\hskip
48.36958pt\times\frac{m-/\\!\\!\\!p}{(m^{2}-p^{2})^{6}}p^{\nu}p^{\sigma}p^{\tau}p^{\lambda}\bigg{]}\bigg{\\}}G_{\rho\sigma;\lambda\tau}G_{\mu\nu},$
(104) $\displaystyle
S_{F}^{6(3)}(0,x)=i\int{\frac{d^{4}p}{(2\pi)^{4}}}e^{-ip\cdot
x}\bigg{(}-\frac{2}{9}\bigg{)}\bigg{\\{}~{}3\gamma^{\sigma}\bigg{[}~{}\gamma^{\rho}\gamma^{\nu}\bigg{(}2\frac{(m+/\\!\\!\\!p){\gamma^{\mu}}}{(m^{2}-p^{2})^{5}}p^{\lambda}p^{\tau}+\frac{{g^{\\{\mu\lambda}}p^{\tau\\}}}{(m^{2}-p^{2})^{4}}\,4p^{\mu}p^{\lambda}p^{\tau}$
$\displaystyle\hskip
48.36958pt\times\frac{1}{(m^{2}-p^{2})^{5}}\bigg{)}+(\mu\leftrightarrow\lambda)+(\rho\leftrightarrow\tau)+(\mu\leftrightarrow\lambda,\rho\leftrightarrow\tau)\bigg{]}+4~{}\gamma^{\rho}\gamma^{\mu}~{}\bigg{[}~{}\frac{m-/\\!\\!\\!p}{(m^{2}-p^{2})^{4}}$
$\displaystyle\hskip 48.36958pt\times
g^{(\nu\lambda\sigma\tau)}~{}+~{}6\frac{m-/\\!\\!\\!p}{(m^{2}-p^{2})^{5}}g^{\\{\nu\lambda}{p^{\sigma}}p^{\tau\\}}~{}+~{}48~{}\frac{m-/\\!\\!\\!p}{(m^{2}-p^{2})^{6}}p^{\nu}p^{\lambda}p^{\sigma}p^{\tau}\bigg{]}~{}\bigg{\\}}~{}G_{\rho\sigma;\tau}G_{\mu\nu;\lambda},$
(105) $\displaystyle
S_{F}^{6(4)}(0,x)=i\int{\frac{d^{4}p}{(2\pi)^{4}}}~{}e^{-ip\cdot
x}~{}\bigg{(}-\frac{1}{2}\bigg{)}\bigg{\\{}\gamma^{\rho}\gamma^{\mu}~{}\bigg{[}~{}\frac{m-/\\!\\!\\!p}{(m^{2}-p^{2})^{4}}g^{\\{\nu\tau\lambda\sigma\\}}+6\frac{m-/\\!\\!\\!p}{(m^{2}-p^{2})^{5}}{g^{\\{\nu\tau}}p^{\lambda}p^{\sigma\\}}$
$\displaystyle\hskip
48.36958pt+48\,\frac{(m-/\\!\\!\\!p)p^{\nu}p^{\tau}p^{\lambda}p^{\sigma}}{(m^{2}-p^{2})^{6}}~{}\bigg{]}+3~{}\gamma^{\sigma}\gamma^{\rho}\gamma^{\nu}~{}\bigg{[}\frac{{g^{\\{\mu\lambda}}p^{\tau\\}}}{(m^{2}-p^{2})^{4}}+\frac{8p^{\mu}p^{\lambda}p^{\tau}}{(m^{2}-p^{2})^{5}}\bigg{]}~{}\bigg{\\}}~{}G_{\rho\sigma}G_{\mu\nu;\lambda\tau},$
(106) $\displaystyle
S_{F}^{6(5)}(0,x)=i\int{\frac{d^{4}p}{(2\pi)^{4}}}e^{-ip\cdot
x}\bigg{\\{}-\frac{i}{18}~{}\gamma^{\nu}\bigg{[}\frac{{(m+/\\!\\!\\!p)}{g^{[\rho\sigma\lambda\tau}}\gamma^{\mu]}}{(m^{2}-p^{2})^{4}}~{}+4~{}\frac{{g^{\\{\rho\sigma\lambda\tau}}p^{\mu\\}}}{(m^{2}-p^{2})^{4}}+6~{}\frac{(m+/\\!\\!\\!p)}{(m^{2}-p^{2})^{5}}$
$\displaystyle\hskip
48.36958pt~{}\cdot~{}g^{\\{\rho\sigma}p^{\lambda}p^{\tau}\gamma^{\mu\\}}+12\frac{{g^{\\{\rho\sigma}}p^{\lambda}p^{\tau}p^{\mu\\}}}{(m^{2}-p^{2})^{5}}-48\frac{{p^{\\{\rho}}p^{\sigma}p^{\lambda}p^{\tau}p^{\mu\\}}(m+/\\!\\!\\!p)}{(m^{2}-p^{2})^{6}}\bigg{]}G_{\mu\nu;\rho\sigma\lambda\tau}\bigg{\\}}.$
(107)

where

$\displaystyle{g^{[\mu\nu\rho\sigma}}p^{\lambda]}=g^{(\mu\nu\rho\sigma)}{p^{\lambda}}+g^{(\lambda\nu\rho\sigma)}p^{\mu}+g^{(\mu\lambda\rho\sigma)}{p^{\nu}}+g^{(\mu\nu\lambda\sigma)}{p^{\rho}}+g^{(\mu\nu\rho\lambda)}{p^{\sigma}},$
(108)
$\displaystyle{g^{\\{\mu\nu}}p^{\rho}p^{\sigma}p^{\lambda\\}}=g^{\mu\nu}p^{\rho}p^{\sigma}p^{\lambda}+g^{\mu\rho}p^{\nu}p^{\sigma}p^{\lambda}+g^{\mu\sigma}p^{\rho}p^{\nu}p^{\lambda}+g^{\mu\lambda}p^{\rho}p^{\sigma}p^{\nu}+g^{\nu\rho}p^{\mu}p^{\sigma}p^{\lambda}$
$\displaystyle\hskip
85.35826pt+g^{\nu\sigma}p^{\rho}p^{\mu}p^{\lambda}+g^{\nu\lambda}p^{\rho}p^{\sigma}p^{\mu}+g^{\rho\sigma}p^{\mu}p^{\nu}p^{\lambda}+g^{\rho\lambda}p^{\mu}p^{\sigma}p^{\nu}+g^{\sigma\lambda}p^{\rho}p^{\mu}p^{\nu},$
(109)
$\displaystyle{g^{\\{\mu\rho}}{\gamma^{\sigma\\}}}=g^{\mu\rho}\gamma^{\sigma}+g^{\mu\sigma}\gamma^{\rho}+g^{\rho\sigma}\gamma^{\mu},$
(110)
$\displaystyle{g^{[\rho\sigma\lambda\tau}}\gamma^{\mu]}=\gamma^{\mu}g^{(\rho\sigma\lambda\tau)}+(\mu\leftrightarrow\rho)+(\mu\leftrightarrow\sigma)+(\mu\leftrightarrow\lambda)+(\mu\leftrightarrow\tau),$
(111)
$\displaystyle{g^{\\{\rho\sigma}}p^{\lambda}\gamma^{\mu\\}}=\gamma^{\mu}g^{\\{\rho\sigma}p^{\lambda\\}}+(\mu\leftrightarrow\rho)+(\mu\leftrightarrow\sigma)+(\mu\leftrightarrow\lambda),$
(112)
$\displaystyle{g^{\\{\rho\sigma}}p^{\lambda}p^{\tau}\gamma^{\mu\\}}=\gamma^{\mu}{g^{\\{\rho\sigma}}p^{\lambda}p^{\tau\\}}+(\mu\leftrightarrow\rho)+(\mu\leftrightarrow\sigma)+(\mu\leftrightarrow\lambda)+(\mu\leftrightarrow\tau),$
(113)
$\displaystyle{\gamma^{\\{\mu}}p^{\rho}p^{\sigma\\}}=\gamma^{\mu}p^{\rho}p^{\sigma}+\gamma^{\rho}p^{\mu}p^{\sigma}+\gamma^{\sigma}p^{\mu}p^{\rho},$
(114)
$\displaystyle{\gamma^{\\{\mu}}p^{\rho}p^{\sigma}p^{\lambda\\}}=\gamma^{\mu}p^{\rho}p^{\sigma}p^{\lambda}+(\mu\leftrightarrow\rho)+(\mu\leftrightarrow\sigma)+(\mu\leftrightarrow\lambda),$
(115)
$\displaystyle{\gamma^{\\{\mu}}p^{\rho}p^{\sigma}p^{\lambda}p^{\tau\\}}=\gamma^{\mu}p^{\rho}p^{\sigma}p^{\lambda}p^{\tau}+(\mu\leftrightarrow\rho)+(\mu\leftrightarrow\sigma)+(\mu\leftrightarrow\lambda)+(\mu\leftrightarrow\tau),$
(116)
$\displaystyle{g^{\mu\\{\lambda}}p^{\tau}\gamma^{\rho\\}}=g^{\mu\lambda}p^{\tau}\gamma^{\rho}+g^{\mu\tau}p^{\lambda}\gamma^{\rho}+g^{\mu\rho}p^{\tau}\gamma^{\lambda}+g^{\mu\tau}p^{\rho}\gamma^{\lambda}+g^{\mu\rho}p^{\lambda}\gamma^{\tau}+g^{\mu\lambda}p^{\rho}\gamma^{\tau}.$
(117)

In these expressions, the mass term is retained explicitly, which application
not only in the case of light quarks, but also in the case of heavy quarks.

Figure 11: Feynman diagrams for the vertex operator
$\Gamma(z\cdot\tensor{D})^{n}$ under the background field theory up to
dimension-six. The symbol “$\times$” attached to the gluon line indicates the
tensor of the local gluon background field, in which “$n$” stands for $n$-th
order covariant derivative.

###  A.2 Vertex operator

The vertex operator $(z\cdot\tensor{D})^{n}$ is encountered when we calculate
the distributed amplitude of moment with the SVZ sum rule. Here it is only the
dimension $d\leq 6$ terms are retained, vertex operator
$(z\cdot\tensor{D})^{n}$ in dimension-six order, the results are as follows

$\displaystyle(z\cdot\tensor{D})_{0}^{n}$
$\displaystyle=(z\cdot\tensor{\partial})^{n},$ (118)
$\displaystyle(z\cdot\tensor{D})_{2}^{n}$
$\displaystyle=-i(z\cdot\tensor{\partial})^{n-1}\underline{x}^{\mu}z^{\nu}G_{\mu\nu},$
(119) $\displaystyle(z\cdot\tensor{D})_{3}^{n}$
$\displaystyle=-\frac{2i}{3}(z\cdot\tensor{\partial})^{n-1}\underline{x}^{\mu}\underline{x}^{\rho}z^{\nu}G_{\mu\nu;\rho},$
(120) $\displaystyle(z\cdot\tensor{D})_{4(1)}^{n}$
$\displaystyle=-\frac{n(n-1)}{2}(z\cdot\tensor{\partial})^{n-2}\underline{x}^{\mu}\underline{x}^{\rho}z^{\nu}z^{\sigma}G_{\mu\nu}G_{\rho\sigma},$
(121) $\displaystyle(z\cdot\tensor{D})_{4(2)}^{n}$
$\displaystyle=\bigg{[}-\frac{i}{4}n(z\cdot\tensor{\partial})^{n-2}\underline{x}^{\mu}\underline{x}^{\rho}\underline{x}^{\sigma}-\frac{i}{12}(n-1)(n-2)(z\cdot\tensor{\partial})^{n-3}\underline{x}^{\mu}z^{\rho}z^{\sigma}z^{\nu}\bigg{]}G_{\mu\nu;\rho\sigma},$
(122) $\displaystyle(z\cdot\tensor{D})_{5(1)}^{n}$
$\displaystyle=-\frac{n(n-1)}{3}(z\cdot\tensor{\partial})^{n-2}\underline{x}^{\mu}\underline{x}^{\rho}\underline{x}^{\sigma}\underline{x}^{\lambda}z^{\nu}z^{\sigma}G_{\mu\nu}G_{\rho\sigma;\lambda},$
(123) $\displaystyle(z\cdot\tensor{D})_{5(2)}^{n}$
$\displaystyle=-\frac{n(n-1)}{3}(z\cdot\tensor{\partial})^{n-2}\underline{x}^{\mu}\underline{x}^{\rho}\underline{x}^{\sigma}\underline{x}^{\lambda}z^{\nu}z^{\sigma}G_{\mu\nu;\lambda}G_{\rho\sigma},$
(124) $\displaystyle(z\cdot\tensor{D})_{5(3)}^{n}$
$\displaystyle=\bigg{[}-\frac{i}{15}n(z\cdot\tensor{\partial})^{n-1}\underline{x}^{\mu}\underline{x}^{\rho}\underline{x}^{\sigma}\underline{x}^{\lambda}z^{\nu}-\frac{i}{45}n(n-1)(n-2)(z\cdot\tensor{\partial})^{n-2}\underline{x}^{\mu}(\underline{x}^{\rho}z^{\sigma}z^{\lambda}$
$\displaystyle+\underline{x}^{\sigma}z^{\rho}z^{\lambda}+\underline{x}^{\lambda}z^{\sigma}z^{\rho})z^{\nu}\bigg{]}G_{\mu\nu;\rho\sigma\lambda},$
(125) $\displaystyle(z\cdot\tensor{D})_{6(1)}^{n}$
$\displaystyle=\frac{i}{6}n(n-1)(n-2)(z\cdot\tensor{\partial})^{n-3}\underline{x}^{\mu}\underline{x}^{\rho}\underline{x}^{\lambda}z^{\nu}z^{\sigma}z^{\tau}G_{\mu\nu}G_{\rho\sigma}G_{\lambda\tau},$
(126) $\displaystyle(z\cdot\tensor{D})_{6(2)}^{n}$
$\displaystyle=\bigg{[}-\frac{n(n-1)}{8}(z\cdot\tensor{\partial})^{n-2}\underline{x}^{\mu}\underline{x}^{\rho}\underline{x}^{\lambda}\underline{x}^{\tau}z^{\nu}z^{\sigma}-\frac{1}{12}n(n-1)(n-2)(n-3)(z\cdot\tensor{\partial})^{n-4}$
$\displaystyle\times\underline{x}^{\mu}\underline{x}^{\rho}z^{\lambda}z^{\tau}z^{\nu}z^{\sigma}\bigg{]}G_{\mu\nu}G_{\rho\sigma;\lambda\tau},$
(127) $\displaystyle\ (z\cdot\tensor{D})_{6(3)}^{n}$
$\displaystyle=-\frac{2}{9}n(n-1)(z\cdot\tensor{\partial})^{n-2}\underline{x}^{\mu}\underline{x}^{\rho}\underline{x}^{\lambda}\underline{x}^{\tau}z^{\nu}z^{\sigma}G_{\mu\nu;\lambda}G_{\rho\sigma;\tau},$
$\displaystyle(z\cdot\tensor{D})_{6(4)}^{n}$
$\displaystyle=-\frac{n(n-1)}{8}(z\cdot\tensor{\partial})^{n-2}\underline{x}^{\mu}\underline{x}^{\rho}\underline{x}^{\lambda}\underline{x}^{\tau}z^{\nu}z^{\sigma}G_{\mu\nu;\lambda\tau}G_{\rho\sigma},$
(128) $\displaystyle(z\cdot\tensor{D})_{6(5)}^{n}$
$\displaystyle=\bigg{[}-\frac{i}{72}\times(z\cdot\tensor{\partial})^{n-1}\underline{x}^{\mu}\underline{x}^{\rho}\underline{x}^{\lambda}\underline{x}^{\tau}z^{\nu}-\frac{i}{216}n(n-1)(n-2)(z\cdot\tensor{\partial})^{n-3}\underline{x}^{\mu}\bigg{(}z^{\rho}z^{\sigma}\underline{x}^{\lambda}\underline{x}^{\tau}$
$\displaystyle+z^{\rho}\underline{x}^{\sigma}z^{\lambda}\underline{x}^{\tau}+z^{\rho}\underline{x}^{\sigma}\underline{x}^{\tau}z^{\tau}+\underline{x}^{\rho}z^{\sigma}z^{\lambda}\underline{x}^{\tau}+\underline{x}^{\rho}z^{\sigma}\underline{x}^{\lambda}z^{\tau}+\underline{x}^{\rho}\underline{x}^{\sigma}z^{\lambda}z^{\tau}\bigg{)}z^{\nu}-\frac{i}{360}n(n-1)$
$\displaystyle\times(n-2)(n-3)(n-4)(z\cdot\tensor{\partial})^{n-5}\underline{x}^{\mu}z^{\rho}z^{\sigma}z^{\lambda}z^{\tau}z^{\nu}\bigg{]}G_{\mu\nu;\lambda\tau\rho\sigma}.$
(129)

where the subscript $k(m)$ with $k=(1,...,6)$ stands for the dimension of the
operator, in which $(m)$ stands for the $m$-th type of the operator with same
dimension. For example, there is two type of dimension four operators, three
type of dimension five operators and five type of dimension-six operators.
Figure 11 represents the Feynman diagrams for the vertex operator
$\Gamma(z\cdot\tensor{D})^{n}$ under the background field theory, where
thirteen figures, i.e. Figs. 11a, $...$, 11m, correspond to
$\Gamma(z\cdot\tensor{D})^{n}_{1},...,\Gamma(z\cdot\tensor{D})^{n}_{6(5)}$,
respectively. The symbol “$\times$” attached to the gluon line indicates the
tensor of the local gluon background field with “$n$” stands for the $n$-th
order covariant derivative.

###  A.3 Vacuum matrix elements

In addition to quark propagators $S_{F}(x,0)$ and vertex operators
$(z\cdot\tensor{D})^{n}$, one may also encounter vacuum condensation in the
process of calculation, which are the basic input parameter. The specific
expression for vacuum matrix element in the $D$-dimension space have the
following forms,

$\displaystyle\langle
0|G^{a}_{\mu\nu}G^{b}_{\rho\sigma}|0\rangle=\frac{1}{96}\delta^{ab}\langle
G^{2}\rangle[g_{\mu\rho}g_{\nu\sigma}-(\sigma\leftrightarrow\rho)],$ (130)
$\displaystyle{\rm Tr}\langle
0|G_{\mu\nu}G_{\rho\sigma}|0\rangle=\frac{\pi}{6}\langle\alpha_{s}G^{2}\rangle[g_{\mu\rho}g_{\nu\sigma}-(\sigma\leftrightarrow\rho)],$
(131) $\displaystyle\langle
0|G^{a}_{\mu\nu;\lambda}G^{b}_{\rho\sigma;\tau}|0\rangle=-\frac{1}{648}\sum_{u,d,s}\delta^{ab}\langle
g_{s}\bar{\psi}\psi\rangle^{2}[2g_{\lambda\tau}(g_{\mu\rho}g_{\nu\sigma}-g_{\mu\sigma}g_{\nu\rho})$
$\displaystyle\hskip
85.35826pt+g_{\lambda\rho}(g_{\tau\mu}g_{\sigma\nu}-g_{\tau\nu}g_{\sigma\mu})-g_{\lambda\sigma}(g_{\tau\mu}g_{\rho\nu}-g_{\tau\nu}g_{\rho\mu})],$
(132) $\displaystyle{\rm Tr}\langle
0|G_{\mu\nu;\lambda}G_{\rho\sigma;\tau}|0\rangle=-\frac{1}{162}\sum_{u,d,s}\langle
g_{s}^{2}\bar{\psi}\psi\rangle^{2}[2g_{\lambda\tau}(g_{\mu\rho}g_{\nu\sigma}-g_{\mu\sigma}g_{\nu\rho})$
$\displaystyle\hskip
85.35826pt+g_{\lambda\rho}(g_{\tau\mu}g_{\sigma\nu}-g_{\tau\nu}g_{\sigma\mu})-g_{\lambda\sigma}(g_{\tau\mu}g_{\rho\nu}-g_{\tau\nu}g_{\rho\mu})],$
(133) $\displaystyle\langle
0|G^{a}_{\mu\nu}G^{b}_{\rho\sigma}G^{c}_{\lambda\tau}|0\rangle=\frac{1}{576}f^{abc}\langle
G^{a}G^{b}G^{c}\rangle\\{[g_{\mu\rho}(g_{\nu\lambda}g_{\sigma\tau}-g_{\nu\tau}g_{\sigma\lambda})$
$\displaystyle\hskip 85.35826pt-
g_{\mu\sigma}(g_{\nu\lambda}g_{\rho\tau}-g_{\nu\tau}g_{\rho\lambda})]-(\mu\leftrightarrow\nu)\\},$
(134) $\displaystyle{\rm Tr}\langle
0|G_{\mu\nu}G_{\rho\sigma}G_{\lambda\tau}|0\rangle=\frac{i}{96}\langle
g_{s}^{3}fG^{3}\rangle\\{[g_{\mu\rho}(g_{\nu\lambda}g_{\sigma\tau}-g_{\nu\tau}g_{\sigma\lambda})$
$\displaystyle\hskip 85.35826pt-
g_{\mu\sigma}(g_{\nu\lambda}g_{\rho\tau}-g_{\nu\tau}g_{\rho\lambda})]-(\mu\leftrightarrow\nu)\\},$
(135) $\displaystyle\langle
0|G^{a}_{\mu\nu}G^{b}_{\rho\sigma;\lambda\tau}|0\rangle=\delta^{ab}\bigg{\\{}\bigg{(}-\frac{1}{1296}\sum_{u,d,s}\langle
g_{s}\bar{\psi}\psi\rangle^{2}-\frac{1}{384}\langle
g_{s}fG^{3}\rangle\bigg{)}$ $\displaystyle\hskip
85.35826pt\times[2g_{\lambda\tau}(g_{\mu\sigma}g_{\nu\rho}-g_{\mu\rho}g_{\nu\sigma})+g_{\rho\tau}(g_{\mu\sigma}g_{\nu\lambda}-g_{\mu\lambda}g_{\nu\sigma})$
$\displaystyle\hskip 85.35826pt-
g_{\sigma\tau}(g_{\mu\lambda}g_{\nu\rho}-g_{\mu\rho}g_{\nu\lambda})]$
$\displaystyle\hskip 85.35826pt+\bigg{(}-\frac{1}{1296}\sum_{u,d,s}\langle
g_{s}\bar{\psi}\psi\rangle^{2}+\frac{1}{384}\langle
g_{s}fG^{3}\rangle\bigg{)}$ $\displaystyle\hskip
85.35826pt\times[g_{\mu\tau}(g_{\rho\nu}g_{\sigma\lambda}-g_{\rho\lambda}g_{\nu\sigma})+(\mu\leftrightarrow\nu)]\bigg{\\}},$
(136) $\displaystyle{\rm Tr}\langle
0|G_{\mu\nu}G_{\rho\sigma;\lambda\tau}|0\rangle=4\bigg{(}-\frac{1}{1296}g_{s}^{2}\sum_{u,d,s}\langle
g_{s}\bar{\psi}\psi\rangle^{2}-\frac{1}{384}\langle
g_{s}^{3}fG^{3}\rangle\bigg{)}$ $\displaystyle\hskip
85.35826pt\times[2g_{\lambda\tau}(g_{\mu\sigma}g_{\nu\rho}-g_{\mu\rho}g_{\nu\sigma})+g_{\rho\tau}(g_{\mu\sigma}g_{\nu\lambda}-g_{\mu\lambda}g_{\nu\sigma})$
$\displaystyle\hskip 85.35826pt-
g_{\sigma\tau}(g_{\mu\lambda}g_{\nu\rho}-g_{\mu\rho}g_{\nu\lambda})]$
$\displaystyle\hskip
85.35826pt+4\bigg{(}-\frac{1}{1296}g_{s}^{2}\sum_{u,d,s}\langle
g_{s}\bar{\psi}\psi\rangle^{2}+\frac{1}{384}\langle
g_{s}^{3}fG^{3}\rangle\bigg{)}$ $\displaystyle\hskip
85.35826pt\times[g_{\mu\tau}(g_{\rho\nu}g_{\sigma\lambda}-g_{\rho\lambda}g_{\nu\sigma})+(\mu\leftrightarrow\nu)],$
(137) $\displaystyle\langle
0|\bar{\psi}_{\alpha}^{a}(x)\psi_{\beta}^{b}(y)|0\rangle=\delta^{ab}\bigg{\\{}\langle\bar{\psi}\psi\rangle\bigg{[}\frac{1}{12}g_{\beta\alpha}+\frac{im_{s}}{48}(/\\!\\!\\!x-/\\!\\!\\!y)_{\beta\alpha}-\frac{m_{s}^{2}}{96}(x-y)^{2}g_{\beta\alpha}$
$\displaystyle\hskip
85.35826pt-\frac{im_{s}^{3}}{576}(x-y)^{2}(/\\!\\!\\!x-/\\!\\!\\!y)_{\beta\alpha}\bigg{]}+g_{s}\langle\bar{\psi}TG\psi\rangle\bigg{[}\frac{1}{192}(x-y)^{2}g_{\beta\alpha}$
$\displaystyle\hskip
85.35826pt+\frac{im_{s}}{864}(x-y)^{2}(/\\!\\!\\!x-/\\!\\!\\!y)_{\beta\alpha}\bigg{]}\bigg{\\}}$
(138) $\displaystyle\langle
0|\bar{\psi}_{\alpha}^{a}(x)\psi_{\beta}^{b}(y)G^{A}_{\mu\nu}|0\rangle=\frac{1}{192}\langle\bar{\psi}TG\psi\rangle(\sigma_{\mu\nu})_{\beta\alpha}(T^{A})^{ba}+\bigg{\\{}-\frac{1}{864}g_{s}\langle\bar{\psi}\psi\rangle^{2}$
$\displaystyle\hskip
85.35826pt\times(g_{\mu\rho}\gamma_{\nu}-g_{\nu\rho}\gamma_{\mu})_{\beta\alpha}(x+y)^{\rho}+i(x-y)^{\rho}\bigg{[}\frac{1}{864}g_{s}\langle\bar{\psi}\psi\rangle^{2}$
$\displaystyle\hskip
85.35826pt+\frac{m_{s}}{384}\langle\bar{\psi}TG\psi\rangle\bigg{]}(\epsilon_{\rho\mu\nu\sigma}\gamma_{5}\gamma^{\sigma})_{\beta\alpha}\bigg{\\}}(T^{A})^{ba}$
(139) $\displaystyle\langle
0|\bar{\psi}_{\alpha}^{a}(0)\psi_{\beta}^{b}(0)G^{A}_{\mu\nu;\rho}|0\rangle=\frac{1}{432}g_{s}\langle\bar{\psi}\psi\rangle^{2}(g_{\mu\rho}\gamma_{\nu}-g_{\nu\rho}\gamma_{\mu})_{\beta\alpha}(T^{A})^{ba}$
(140)

the symbol “${\rm Tr}$” stands for tracing to the colour matrixes. Obviously,

$\displaystyle\langle
0|G^{B}_{\rho\sigma;\lambda\tau}G^{A}_{\mu\nu}|0\rangle=\langle
0|G^{A}_{\mu\nu}G^{B}_{\rho\sigma;\lambda\tau}|0\rangle.$ (141)

As a final remark, by making use of the equation

$[\widetilde{D}_{\mu},\widetilde{D}_{\nu}]^{AB}=-g_{s}f^{ABC}G^{C}_{\mu\nu},$
(142)

we can obtain a useful relation

$\langle g_{s}^{3}fG^{3}\rangle=\frac{8}{27}\sum_{u,d,s}\langle
g_{s}^{2}\bar{\psi}\psi\rangle^{2}.$ (143)

which, can also be referred to our previous work Zhong:2014jla .

##  Appendix B Detailed expressions for $I_{m_{s}^{0}}(n,M^{2})$ and
$I_{m_{s}^{2}}(n,M^{2})$ in BFTSR for
$\langle\xi_{2;\eta^{(\prime)}}^{n}\rangle|_{\mu}$

According to the Feynman rule, we can obtain the following detailed
expressions within $I_{m_{s}^{0}}(n,M^{2})$ and $I_{m_{s}^{2}}(n,M^{2})$ for
the final BFTSR, i.e. Eq. (39). Some of the terms have the symmetric
structure, we can combined these terms together. Here, we only list the key
expression for the last step of the BFTSR. Firstly, we give the
$I_{m_{s}^{0}}(n,M^{2})$ terms which are expressed as,

$\displaystyle\hat{\cal
B}_{M^{2}}I_{000}=\frac{3}{8\pi^{2}}\frac{1+(-1)^{n}}{(n+1)(n+3)},$ (144)
$\displaystyle\hat{\cal B}_{M^{2}}I_{006+600}=~{}\frac{g_{s}^{2}\sum\langle
g_{s}\bar{q}q\rangle^{2}}{M^{6}}~{}[1+(-1)^{n}]~{}~{}\frac{-17}{324\pi^{2}}~{}\bigg{\\{}~{}\left[(2n+1)~{}\left(-\ln\frac{M^{2}}{\mu^{2}}\right)~{}-(n+2)\right]$
$\displaystyle\hskip
59.75095pt+\theta(n-2)~{}\bigg{[}~{}\frac{1}{2}~{}(2n+1)~{}\bigg{(}\psi\Big{(}\frac{n+1}{2}\Big{)}~{}-~{}\psi\Big{(}\frac{n}{2}\Big{)}~{}+~{}\ln
4\bigg{)}~{}-\frac{(n+1)^{2}}{n}\bigg{]}~{}\bigg{\\}},$ (145)
$\displaystyle\hat{\cal B}_{M^{2}}I_{024+420}=\frac{g_{s}^{2}\sum\langle
g_{s}\bar{q}q\rangle^{2}}{M^{6}}[1+(-1)^{n}]\frac{-1}{162\pi^{2}}\theta(n-2)\bigg{\\{}n\bigg{(}-\ln\frac{M^{2}}{\mu^{2}}\bigg{)}+\bigg{[}\frac{n}{2}\bigg{(}\psi\Big{(}\frac{n+1}{2}\Big{)}$
$\displaystyle\hskip 59.75095pt-\psi\Big{(}\frac{n}{2}\Big{)}+\ln
4\bigg{)}-1\bigg{]}\bigg{\\}},$ (146) $\displaystyle\hat{\cal
B}_{M^{2}}I_{033+330}=\frac{g_{s}^{2}\sum\langle
g_{s}\bar{q}q\rangle^{2}}{M^{6}}~{}[1+(-1)^{n}]~{}\frac{2}{243\pi^{2}}~{}\theta(n-2)~{}\bigg{\\{}~{}n~{}\left(-\ln\frac{M^{2}}{\mu^{2}}~{}-~{}\frac{1}{2}\right)~{}+\bigg{[}\frac{n}{2}$
$\displaystyle\hskip
59.75095pt\times\left(\psi\Big{(}\frac{n+1}{2}\Big{)}-\psi\Big{(}\frac{n}{2}\Big{)}+\ln
4\right)-1\bigg{]}\bigg{\\}},$ (147) $\displaystyle\hat{\cal
B}_{M^{2}}I_{202+040}=\frac{\langle\alpha_{s}G^{2}\rangle}{M^{4}}\frac{1+n\theta(n-2)}{24\pi(n+1)}[1+(-1)^{n}],$
(148) $\displaystyle\hat{\cal B}_{M^{2}}I_{060}=\frac{16g_{s}^{2}\sum\langle
g_{s}\bar{q}q\rangle^{2}-81\langle
g_{s}^{3}fG^{3}\rangle}{M^{6}}\frac{n}{7776\pi^{2}}[1+(-1)^{n}]\,\theta(n-2),$
(149) $\displaystyle\hat{\cal B}_{M^{2}}I_{204+402}=\frac{g_{s}^{2}\sum\langle
g_{s}\bar{q}q\rangle^{2}}{M^{6}}[1+(-1)^{n}]\,\frac{1}{324\pi^{2}}\bigg{\\{}\left(-\ln\frac{M^{2}}{\mu^{2}}\right)+\theta(n-2)~{}\bigg{[}~{}\frac{1}{2}\bigg{(}\psi\Big{(}\frac{n+1}{2}\Big{)}$
$\displaystyle\hskip 59.75095pt-\psi\Big{(}\frac{n}{2}\Big{)}+\ln
4\bigg{)}-\frac{1}{n}\bigg{]}\bigg{\\}},$ (150) $\displaystyle\hat{\cal
B}_{M^{2}}I_{303(1)+303(3)}=~{}\frac{g_{s}^{2}\sum\langle
g_{s}\bar{q}q\rangle^{2}}{M^{6}}~{}[1+(-1)^{n}]~{}~{}\frac{-2}{243\pi^{2}}~{}\bigg{\\{}~{}\bigg{(}~{}-\ln\frac{M^{2}}{\mu^{2}}~{}\bigg{)}~{}+\theta(n-2)~{}\bigg{[}~{}\frac{1}{2}$
$\displaystyle\hskip
79.6678pt\times\bigg{(}\psi\Big{(}\frac{n+1}{2}\Big{)}-\psi\Big{(}\frac{n}{2}\Big{)}+\ln
4\bigg{)}-\frac{1}{n}\bigg{]}\bigg{\\}},$ (151) $\displaystyle\hat{\cal
B}_{M^{2}}I_{303(2)}=\frac{g_{s}^{2}\sum\langle
g_{s}\bar{q}q\rangle^{2}}{M^{6}}[1+(-1)^{n}]\frac{1}{162\pi^{2}}~{}\bigg{\\{}\bigg{(}-\ln\frac{M^{2}}{\mu^{2}}+\frac{1}{2}\bigg{)}+\frac{\theta(n-2)}{2}\bigg{(}\,\psi\Big{(}\frac{n+1}{2}\Big{)}$
$\displaystyle\hskip 51.21504pt-\psi\Big{(}\frac{n}{2}\Big{)}+\ln
4\bigg{)}\bigg{\\}}.$ (152)

In which, the subscript $a,b,c,m$ in $\hat{\cal B}_{M^{2}}I_{abc(m)}$ for Eqs.
(144)-(152) have the meaning that, $a,b,c$ stand for $a,b,c$th-order of the
propagator, vertex and propagator, while $m$ means the $m$th terms of
$I^{m_{s}^{2}}_{abc}$. Then, we give the expression for
$I_{m_{s}^{2}}(n,M^{2})$ terms separately,

$\displaystyle\hat{\cal
B}_{M^{2}}I_{004+400}^{m_{s}^{2}}=\frac{\langle\alpha_{s}G^{2}\rangle}{M^{6}}m_{s}^{2}\frac{-1}{12\pi}[1+(-1)^{n}]\bigg{\\{}~{}\bigg{[}~{}(2n+1)\Big{(}-\ln\frac{M^{2}}{\mu^{2}}+\psi(n+1)+2\gamma_{E}\Big{)}-n$
$\displaystyle\hskip
59.75095pt-2~{}\bigg{]}~{}+~{}\theta(n-2)\bigg{[}~{}\frac{2n+1}{2}\bigg{(}\psi\Big{(}\frac{n+1}{2}\Big{)}-\psi\Big{(}\frac{n}{2}\Big{)}+\ln
4\bigg{)}-\frac{(n+1)^{2}}{n}\bigg{]}\bigg{\\}},$ (153)
$\displaystyle\hat{\cal B}_{M^{2}}I_{006(1)+600(1)}^{m_{s}^{2}}=\frac{\langle
g_{s}^{2}\bar{q}q\rangle^{2}}{M^{8}}m_{s}^{2}\frac{17}{324\pi^{2}}[1+(-1)^{n}]\bigg{\\{}\delta^{n0}+2n\Big{(}-\ln\frac{M^{2}}{\mu^{2}}+\psi(n+1)+2\gamma_{E}\Big{)}$
$\displaystyle\hskip
79.6678pt-n-3+\theta(n-2)~{}\left[~{}n\left(\psi\Big{(}\frac{n+1}{2}\Big{)}~{}-~{}\psi\Big{(}\frac{n}{2}\Big{)}~{}+~{}\ln
4\right)-\left(n+1\right)\,\right]~{}\bigg{\\}},$ (154)
$\displaystyle\hat{\cal B}_{M^{2}}I_{006(2)+600(2)}^{m_{s}^{2}}=\frac{8\langle
g_{s}^{2}\bar{q}q\rangle^{2}-3\langle
g_{s}^{3}fG^{3}\rangle}{M^{8}}m_{s}^{2}\frac{1}{144\pi^{2}}[1+(-1)^{n}]\,\bigg{\\{}\frac{1}{2}~{}\delta^{n0}-\bigg{[}\,2n^{2}\bigg{(}-\ln\frac{M^{2}}{\mu^{2}}$
$\displaystyle\hskip
79.6678pt+\psi(n+1)+2\gamma_{E}\bigg{)}-\frac{7n^{2}+19n+2}{4}\bigg{]}\bigg{\\}}+\bigg{\\{}\theta(n-3)\bigg{[}-n^{2}\bigg{(}\,\psi\Big{(}\frac{n+1}{2}\Big{)}$
$\displaystyle\hskip
79.6678pt-\psi\Big{(}\frac{n}{2}\Big{)}+\frac{1}{4}(5n^{2}+5n+2)\bigg{]}\bigg{\\}},$
(155) $\displaystyle\hat{\cal
B}_{M^{2}}I_{024(1)+420(1)}^{m_{s}^{2}}=\frac{\langle
g_{s}^{2}\bar{q}q\rangle^{2}}{M^{8}}m_{s}^{2}\frac{n}{162\pi^{2}}[1+(-1)^{n}]\,\bigg{\\{}\theta(n-1)\bigg{(}-\ln\frac{M^{2}}{\mu^{2}}+\psi(n+1)+2\gamma_{E}$
$\displaystyle\hskip
79.6678pt-\frac{1}{n}\bigg{)}+\frac{1}{2}\theta(n-2)\left[\psi\Big{(}\frac{n+1}{2}\Big{)}-\psi\Big{(}\frac{n}{2}\Big{)}+\ln
4\right]\bigg{\\}},$ (156) $\displaystyle\hat{\cal
B}_{M^{2}}I_{024(2)+420(2)}^{m_{s}^{2}}=\frac{8\langle
g_{s}^{2}\bar{q}q\rangle^{2}-27\langle
g_{s}^{3}fG^{3}\rangle}{M^{8}}~{}m_{s}^{2}~{}\frac{-n}{3888\pi^{2}}~{}[1+(-1)^{n}]\,~{}\bigg{\\{}~{}\theta(n-1)~{}\bigg{[}(2n-1)$
$\displaystyle\hskip
79.6678pt\times\bigg{(}-\ln\frac{M^{2}}{\mu^{2}}+\psi(n+1)+2\gamma_{E}~{}\bigg{)}-n-3~{}+\frac{1}{n}~{}\bigg{]}+~{}\theta(n-3)~{}\bigg{[}~{}\frac{2n-1}{2}$