The ratio $\mathcal{R}(D_{s})$ for $B_{s}\to D_{s}\ell\nu_{\ell}$ by using the QCD light-cone sum rules within the framework of heavy quark effective field theory

For the $B_{s}\to D_{s}l\bar{\nu}_{l}$ decays, the transition matrix element can be parameterized as follows:

			$\displaystyle\langle D_{s}(p)|\bar{c}\gamma_{\mu}b|B_{s}(p+q)\rangle$		(1)
		$\displaystyle=$	$\displaystyle 2f_{+}^{B_{s}\to D_{s}}(q^{2})p_{\mu}+\left[f_{+}^{B_{s}\to D_{s}}(q^{2})+f_{-}^{B_{s}\to D_{s}}(q^{2})\right]q_{\mu}$

and

$\displaystyle f_{0}^{B_{s}\to D_{s}}(q^{2})=f_{+}^{B_{s}\to D_{s}}(q^{2})+\frac{q^{2}}{m_{B_{s}}^{2}-m_{D_{s}}^{2}}f_{-}^{B_{s}\to D_{s}}(q^{2}),$

(2)

where $p$ is the momentum of the $D_{s}$-meson and $(p+q)$ is the momentum of $B_{s}$-meson. At the maximum recoil point, we have $f_{+}^{B_{s}\to D_{s}}(0)=f_{0}^{B_{s}\to D_{s}}(0)$. The transition matrix element can be expanded as $1/m_{b}$-power series within the framework of HQEFT. Based on the heavy quark symmetry, the transition matrix element of heavy quark in the effective theory is parameterized as Wang:1999zd ; Wang:2000sc ; Wang:2000gs :

	$\displaystyle\langle D_{s}(p)|\bar{c}\gamma_{\mu}b|B_{s}(p+q)\rangle$	$\displaystyle=$	$\displaystyle\frac{\sqrt{m_{B_{s}}}}{\sqrt{\bar{\Lambda}_{B_{s}}}}\langle D_{s}(p)|\bar{u}\gamma_{\mu}b_{v}^{+}|{B_{s}}_{v}\rangle$		(3)
		$\displaystyle=$	$\displaystyle-\frac{\sqrt{m_{B_{s}}}}{\sqrt{\bar{\Lambda}_{B_{s}}}}{\rm Tr}[D_{s}(v,p){\gamma_{\mu}}{\mathcal{M}}_{v}],$

where

	$\displaystyle\bar{\Lambda}_{B_{s}}$	$\displaystyle=$	$\displaystyle m_{B_{s}}-m_{b},$
	$\displaystyle D_{s}(v,p)$	$\displaystyle=$	$\displaystyle{\gamma^{5}}[A(v\cdot p)+\not\!\hat{p}B(v\cdot p)],$		(4)
	$\displaystyle{\mathcal{M}}_{v}$	$\displaystyle=$	$\displaystyle-\sqrt{\bar{\Lambda}}(1+\not\!v)\gamma^{5}/2,$

where $b_{v}^{+}$ is the effective $b$-quark field and $v$ is the $B_{s}$-meson velocity, $\hat{p}^{\mu}=p^{\mu}/(v\cdot p)$. $A(v\cdot p)$ and $B(v\cdot p)$ are leading-order heavy flavor-spin independent coefficient functions. $\bar{\Lambda}=\mathop{\lim}\limits_{m_{b}\to\infty}{\bar{\Lambda}_{B_{s}}}$, which is the heavy flavor independent binding energy that reflects the effects of light degrees of freedom in the heavy hadron. Using those formulas, we obtain the $B_{s}\to D_{s}\ell\bar{\nu}_{\ell}$ TFFs $f_{\pm}(q^{2})$, which are

$\displaystyle f_{\pm}^{B_{s}\to D_{s}}(q^{2})=\frac{\sqrt{\bar{\Lambda}}}{\sqrt{m_{B_{s}}\bar{\Lambda}_{B_{s}}}}\left[A(y)\pm\frac{m_{B_{s}}}{y}B(y)\right]+\cdots,$

(5)

with

$$y=v\cdot p=(m^{2}_{B_{s}}+m^{2}_{D_{s}}-q^{2})/(2m_{B_{s}}),$$

(6)

where “$\cdots$” denotes the higher-order ${\cal O}(1/{m_{b}})$ contributions that will not be taken into consideration here.

To derive the sum rules of the two leading order heavy flavor-spin independent coefficient functions $A(y)$ and $B(y)$, we construct the following correlator:

$\displaystyle F_{\mu}(p,q)=i\int d^{4}xe^{iq\cdot x}\langle D_{s}(p)|T\{j_{n}(x),j^{\dagger}_{0}(0)\}|0\rangle,$

(7)

where the currents

	$\displaystyle j_{n}(x)$	$\displaystyle=$	$\displaystyle\bar{c}(x)\gamma_{\mu}(1+\gamma_{5})b(x)$		(8)
	$\displaystyle j^{{\dagger}}_{0}(0)$	$\displaystyle=$	$\displaystyle\bar{b}(0)i(1+\gamma_{5})s(0)$		(9)

Following the standard procedure of LCSR approach, we first deal with the hadronic representation for the correlation function. One can insert a complete series of the intermediate hadronic states in the correlator (15) in the physical $q^{2}$-region and isolate the pole term of the lowest pseudoscalar state from the hadronic representation. Then the correlator $F_{\mu}(p,q)$ becomes:

	$\displaystyle F_{\mu}^{\rm Had.}(p,q)=\frac{\langle D_{s}(p)|\bar{c}\gamma_{\mu}b|B_{s}\rangle\langle B_{s}|\bar{b}i\gamma_{5}s|0\rangle}{m_{B_{s}}^{2}-(p+q)^{2}}$		(10)
	$\displaystyle\quad+\sum\limits_{B_{s}^{H}}\frac{\langle D_{s}(p)|\bar{c}\gamma_{\mu}(1+\gamma_{5})b|B_{s}^{H}\rangle\langle B_{s}^{H}|\bar{b}i(1+\gamma_{5})s|0\rangle}{m_{B_{s}^{H}}^{2}-(p+q)^{2}}.$
			(11)

In the effective theory of heavy quark, the hadronic representation (11) can be expanded in powers of $1/m_{b}$. Taking the transition matrix element (3) into consideration and neglecting the contributions from higher $1/m_{b}$ order, we can further write the hadronic representation as:

	$\displaystyle F_{\mu}^{\rm Had.}(p,q)$	$\displaystyle=$	$\displaystyle 2F\frac{A(y)v^{\mu}+B(y)\hat{p}^{\mu}}{2\bar{\Lambda}_{B_{s}}-2v\cdot k}$		(13)
			$\displaystyle+\int_{s_{0}}^{\infty}ds\frac{\rho(y,s)}{s-2v\cdot k}+{\rm Subtractions},$

with the matrix element Korner:2002ba

$\displaystyle\langle B_{s}|\bar{b}^{+}_{v}i\gamma_{5}d|0\rangle=\frac{i}{2}F{\rm Tr}[\gamma_{5}\mathcal{M}_{v}],$

(14)

where $F$ is the leading-order decay constant of the $B_{s}$-meson Wang:2001mi ; Wang:2002zba . $k$ is the residual momentum of the heavy hadronic. Using the ansatz of the quark-hadron duality the spectral density $\rho(y,s)$ can be obtained Shifman:1978bx ; Shifman:1978by .

On the other hand, we apply the operator product expansion (OPE) to the correlator in the deep Euclidean region. The correlator (15) can be explicitly written as

	$\displaystyle F_{\mu}(p,q)$	$\displaystyle=$	$\displaystyle i\int d^{4}xe^{i(q-m_{b})\cdot x}\langle D_{s}(p)|T\{\bar{c}(x)\gamma_{\mu}(1+\gamma_{5})b^{+}_{v}(x),$		(15)
			$\displaystyle\bar{b}^{+}_{v}(0)i(1+\gamma_{5})s(0)\}|0\rangle.$

Using the $B$-meson heavy-quark propagator $S(x,v)=(1+/\!\!\!v)\times\int_{0}^{\infty}dt\delta(x-vt)/2$ Wang:2001mi , the correlator can be expanded as a complex power series over the $D_{s}$-meson LCDAs. Due to the chiral suppressions, it is noted that the main contribution to the correlator comes from the leading-twist LCDA, and the contributions from all the twist-3 LCDAs are exactly zero.

Through the dispersion relation, the OPE in deep Euclidean region and the hadron expression in physical region can be matched. And by further applying the Borel transformation to suppress the contributions from power-suppressed terems, the LCSRs for the coefficient functions $A(y)$ and $B(y)$ are

	$\displaystyle A(y)$	$\displaystyle=$	$\displaystyle-\frac{f_{D_{s}}}{2F}\int_{0}^{s_{0}^{B_{s}}}dse^{(2\bar{\Lambda}_{B_{s}}-s)/T}\frac{1}{y^{2}}\frac{\partial}{\partial u}g_{2}(u)\bigg{|}_{u=1-\frac{s}{2y}},$		(16)
	$\displaystyle B(y)$	$\displaystyle=$	$\displaystyle-\frac{f_{D_{s}}}{2F}\int_{0}^{s_{0}^{B_{s}}}dse^{(2\bar{\Lambda}_{B_{s}}-s)/T}\bigg{[}-\phi_{2;D_{s}}(u)$		(17)
			$\displaystyle+\bigg{(}\frac{1}{y}\frac{\partial}{\partial u}\bigg{)}^{2}g_{1}(u)-\frac{1}{y^{2}}\frac{\partial}{\partial u}g_{2}(u)\bigg{]}\bigg{|}_{u=1-\frac{s}{2y}},$

where $T$ is the Borel parameter and $s_{0}^{B_{s}}$ is the continuum threshold, $g_{1}$ and $g_{2}$ are twist-four LCDAs. Since the contributions from the twist-4 LCDAs are only several percent, so we shall directly adopt the light pseudo-scale ones to do the calculations, whose explicit forms can be found in Refs.Ball:1998je . Substituting them into Eqs.(2,5), we obtain

	$\displaystyle f_{+}^{B_{s}\to D_{s}}(q^{2})$	$\displaystyle=$	$\displaystyle-\frac{f_{D_{s}}\sqrt{\bar{\Lambda}}}{2F\sqrt{m_{B_{s}}\bar{\Lambda}_{B_{s}}}}\int_{0}^{s_{0}^{B_{s}}}dse^{(2\bar{\Lambda}_{B_{s}}-s)/T}$		(18)
		$\displaystyle\times$	$\displaystyle\bigg{\{}\frac{1}{y^{2}}\frac{\partial}{\partial u}g_{2}(u)+\frac{m_{B_{s}}}{y}\bigg{[}-\phi_{2;D_{s}}(u)$
		$\displaystyle+$	$\displaystyle\bigg{(}\frac{1}{y}\frac{\partial}{\partial u}\bigg{)}^{2}g_{1}(u)-\frac{1}{y^{2}}\frac{\partial}{\partial u}g_{2}(u)\bigg{]}\bigg{\}}\bigg{|}_{u=1-\frac{s}{2y}},$
	$\displaystyle f_{0}^{B_{s}\to D_{s}}(q^{2})$	$\displaystyle=$	$\displaystyle-\frac{f_{D_{s}}\sqrt{\bar{\Lambda}}}{2F\sqrt{m_{B_{s}}\bar{\Lambda}_{B_{s}}}}\int_{0}^{s_{0}^{B_{s}}}dse^{(2\bar{\Lambda}_{B_{s}}-s)/T}$		(19)
		$\displaystyle\times$	$\displaystyle\bigg{\{}\bigg{(}1+\frac{q^{2}}{m_{B_{s}}^{2}-m_{D_{s}}^{2}}\bigg{)}\frac{1}{y^{2}}\frac{\partial}{\partial u}g_{2}(u)$
		$\displaystyle+$	$\displaystyle\bigg{(}1-\frac{q^{2}}{m_{B_{s}}^{2}-m_{D_{s}}^{2}}\bigg{)}\frac{m_{B_{s}}}{y}\bigg{[}-\phi_{2;D_{s}}(u)$
		$\displaystyle+$	$\displaystyle\bigg{(}\frac{1}{y}\frac{\partial}{\partial u}\bigg{)}^{2}g_{1}(u)-\frac{1}{y^{2}}\frac{\partial}{\partial u}g_{2}(u)\bigg{]}\bigg{\}}\bigg{|}_{u=1-\frac{s}{2y}}.$

The Borel parameter $T$ and the continuum threshold $s_{0}^{B_{s}}$ shall be fixed such that the resulting TFFs do not depend too much on the precise values of those parameters. In addition, the continuum contribution, which is the part of the dispersive integral from $s_{0}^{B_{s}}$ to $\infty$ that is subtracted from both sides of the equation, should not be too large.

To determine the TFFs $f_{+,0}^{B_{s}\to D_{s}}(q^{2})$ of the exclusive process $B_{s}\to D_{s}\ell\bar{\nu}_{\ell}$, we take Zyla:2020zbs ; Wang:2000sc

	$\displaystyle m_{B_{s}}$	$\displaystyle=$	$\displaystyle 5.367\pm 0.00014{\rm GeV},$
	$\displaystyle m_{D_{s}}$	$\displaystyle=$	$\displaystyle 1.968\pm 0.00007{\rm GeV},$
	$\displaystyle f_{D_{s}}$	$\displaystyle=$	$\displaystyle 0.256\pm 0.0042{\rm GeV},$
	$\displaystyle F$	$\displaystyle=$	$\displaystyle 0.30\pm 0.04{\rm GeV^{3/2}}.$

For the $D_{s}$-meson leading-twist LCDA $\phi_{2;D_{s}}(u,\mu)$ in the LCSRs of TFFs (18) and (19), we adopt the light cone harmonic oscillator (LCHO) model suggested in Ref.Zhang:2021wnv , which has a better end-point behavior and takes the following form:

	$\displaystyle\phi_{2;D_{s}}(u,\mu)$	$\displaystyle=$	$\displaystyle\frac{\sqrt{6}A_{D_{s}}\beta_{D_{s}}^{2}}{\pi^{2}f_{D_{s}}}u\bar{u}\varphi_{2;D_{s}}(u)$		(20)
		$\displaystyle\times$	$\displaystyle\exp\left[-\frac{\hat{m}_{c}^{2}u+\hat{m}_{s}^{2}\bar{u}}{8\beta_{D_{s}}^{2}u\bar{u}}\right]$
		$\displaystyle\times$	$\displaystyle\left\{1-\exp\left[-\frac{\mu^{2}}{8\beta_{D_{s}}^{2}u\bar{u}}\right]\right\},$

where $\bar{u}=1-u$ and $\varphi_{2;D_{s}}(u)=1+\sum^{4}_{n=1}B_{n}^{D_{s}}\times C_{n}^{3/2}(\xi)$ with $\xi=u-\bar{u}$. $\hat{m}_{c}$ and $\hat{m}_{s}$ are $c$- and $s$-constituent quark masses, whose values are taken as $\hat{m}_{c}\simeq 1.5~{}{\rm GeV}$ and $\hat{m}_{s}\simeq 0.5~{}{\rm GeV}$. $A_{D_{s}}$. $\beta_{D_{s}}$, $B_{1}^{D_{s}}$, $B_{2}^{D_{s}}$, $B_{3}^{D_{s}}$ and $B_{4}^{D_{s}}$ are model parameters, whose initial values at the scale $\mu=2{\rm GeV}$ have been given in Ref. Zhang:2021wnv . For the present process, the typical factorization scale $\mu$ = $(m^{2}_{B_{s}}-m_{b}^{2})^{1/2}=3{\rm GeV}$. The input parameters at the scale $\mu=3{\rm GeV}$ can be achieved by using the conventional one-loop evolution equation Lepage:1980fj , and these values are given in Table 1. Figure 1 shows the behavior of the $D_{s}$-meson leading-twist LCDA $\phi_{2;D_{s}}(u,\mu=3{\rm GeV})$ with the typical values exhibited in Table 1, where the solid line is the central value and the shaded band shows its uncertainty given in Table 1.

Next, we calculate the $B_{s}\to D_{s}$ TFFs $f^{B_{s}\to D_{s}}_{+}(q^{2})$ and $f^{B_{s}\to D_{s}}_{0}(q^{2})$ by using the LCSRs (18) and (19) in the $q^{2}$-region when the LCSR approach is applicable, i.e., $0<q^{2}<7{\rm GeV^{2}}$. For the purpose, we first determine the continuum threshold $s_{0}^{B_{s}}$ and the Borel parameter $T$. Within the framework of HQEFT, the continuum threshold $s_{0}^{B_{s}}\equiv 2\bar{\Lambda}_{B_{s}^{*}}=2(m_{B_{s}^{*}}-m_{b})$ with $B_{s}^{*}$ being the $B_{s}$-meson first excited state Wang:2001mi , and we take $s_{0}^{B_{s}}=3.85\pm 0.15{\rm GeV}$. To determine the Borel window, as suggested in Refs.Wang:2001mi ; Wang:2002zba , we require the TFFs $f^{B_{s}\to D_{s}}_{+,0}(q^{2})$ to be as stable as possible within corresponding Borel windows. Figure 2 shows the TFFs $f^{B_{s}\to D_{s}}_{+,0}(q^{2})$ versus the Borel parameter $T$ at several typical squared momenta transfer, in which the solid, the dashed, the dot-dashed and the dotted lines are for $q^{2}=0,3,5,7{\rm GeV}^{2}$, respectively. One can find that, the TFFs $f^{B_{s}\to D_{s}}_{+,0}(q^{2})$ are stable for a large $T$, e.g. the uncertainty caused by $T$ is less than $5\%$ when $T\geq 10{\rm GeV}$ for all those $q^{2}$ values. Therefore, we take the Borel window as $10{\rm GeV}<T<20{\rm GeV}$.

At the maximum recoil point $q^{2}=0$, we have

	$\displaystyle f^{B_{s}\to D_{s}}_{+,0}(0)$	$\displaystyle=$	$\displaystyle 0.533^{+0.082}_{-0.063}|_{\phi_{2;D_{s}}}\ ^{+0.007}{}_{-0.014}|_{T}\ ^{+0.065}{}_{-0.063}|_{s_{0}^{B_{s}}}$		(21)
			${}^{+0.027}_{-0.037}|_{F}\ ^{+0.004}{}_{-0.004}|_{m_{b}}.$

There are also errors caused by the uncertainties of $m_{B_{s}}$ and $m_{D_{s}}$, which are negligibly small. It is fond that the uncertainties of the $D_{s}$-meson leading-twist LCDA $\phi_{2;D_{s}}$ and the continuum threshold $s_{0}^{B_{s}}$ are main errors of $f^{B_{s}\to D_{s}}_{+,0}(0)$. By adding all those errors in quadrature, we obtain $f^{B_{s}\to D_{s}}_{+,0}(0)=0.533^{+0.160}_{-0.128}$. We present the theoretical predictions of $f_{+,0}^{B_{s}\to D_{s}}(q^{2})$ at the maximum recoil point $q^{2}=0$ in Table 2, where the predictions under the pQCD approach Fan:2013kqa , the pQCD+LQCD approach Hu:2019bdf , the LQCD approach Monahan:2017uby ; Monahan:2018lzv , the QCD SR approach Blasi:1993fi ; Azizi:2008tt , the LCSR approach Li:2009wq , the BSE approach Chen:2011ut and the RQM approach Bhol:2014jta are also presented. Our present prediction of $f^{B_{s}\to D_{s}}_{+,0}(0)$ is in good agreement with the values calculated with the pQCD prediction Fan:2013kqa , the pQCD+LQCD prediction Hu:2019bdf and the BSE prediction Chen:2011ut .

As mentioned above, the LCSRs (18) and (19) for TFFs $f^{B_{s}\to D_{s}}_{+,0}(q^{2})$ are only reliable in low and intermediate regions, i.e., $0<q^{2}<7{\rm GeV^{2}}$. To estimate the total decay width of the semi-leptonic decay $B_{s}\to D_{s}\ell\bar{\nu}_{\ell}$, we extrapolate the TFFs to the whole physically allowable $q^{2}$-region, $0<q^{2}<(m_{B_{s}}-m_{D_{s}})^{2}=11.50{\rm GeV^{2}}$, via the $z$-series parametrization Khodjamirian:2011ub ; Bourrely:2008za :

	$\displaystyle f_{+}^{B_{s}\to D_{s}}(q^{2})$	$\displaystyle=$	$\displaystyle\frac{f_{+}^{B_{s}\to D_{s}}(0)}{1-q^{2}/m_{B_{s}^{*}}^{2}}\bigg{\{}1+\sum_{k=1}^{N-1}b_{k}\bigg{[}z(q^{2})^{k}-z(0)^{k}$		(22)
			$\displaystyle-(-1)^{N-k}\frac{k}{N}(z(q^{2})^{N}-z(0)^{N})\bigg{]}\bigg{\}},$
	$\displaystyle f_{0}^{B_{s}\to D_{s}}(q^{2})$	$\displaystyle=$	$\displaystyle f_{0}^{B_{s}\to D_{s}}(0)\bigg{\{}1+\sum_{k=1}^{N-1}b_{k}(z(q^{2})^{k}-z(0)^{k})\bigg{\}}.$

where

	$\displaystyle z(q^{2})$	$\displaystyle=$	$\displaystyle\frac{\sqrt{t_{+}-q^{2}}-\sqrt{t_{+}-t_{0}}}{\sqrt{t_{+}-q^{2}}+\sqrt{t_{+}-t_{0}}},$
	$\displaystyle t_{0}$	$\displaystyle=$	$\displaystyle t_{+}(1-\sqrt{1-t_{+}/t_{-}}),$
	$\displaystyle t_{\pm}$	$\displaystyle=$	$\displaystyle(m_{B_{s}}\pm m_{D_{s}})^{2}.$

Then, by fitting the values of the TFFs in low and intermediate regions calculated via the LCSRs (18) and (19), the coefficients $b_{1}$, $b_{2}$ and $b_{3}$ in extrapolation formula (22) and (LABEL:Eq:f0EX) can be determined, and which have been exhibited in Table 3. The quality-of-fit is defined as:

$$\Delta=\frac{\sum_{t}{|f^{i}-f^{\rm fit}|}}{\sum_{t}{|f^{i}|}}\times 100\%,t\in\bigg{\{}0,\frac{1}{2},...,\frac{23}{2},12\bigg{\}}~{}\rm{GeV}^{2}.$$

(24)

The coefficients $b_{i}$ are determined such that the quality-of-fit $(\Delta)$ is no more than $1\%$. The $\Delta$ values for the central, the upper and lower TFFs are shown in Table 3. These quality-of-fits are much smaller than $1\%$, indicating that our present extrapolations are of high accuracy. We present the extrapolated TFFs $f^{B_{s}\to D_{s}}_{+,0}(q^{2})$ in Figure 3, where the shaded hands are theoretical uncertainties from all the mentioned error sources. For comparison, we present the results of the pQCD+LQCD approach Hu:2019bdf , the pQCD approach Hu:2019bdf , the RQM approach Faustov:2012mt and the LQCD approach Monahan:2017uby .

The branching fraction of the semileptonic decay $B_{s}\to D_{s}\ell\bar{\nu}_{\ell}$ is defined as

$$\mathcal{B}(B_{s}\to D_{s}\ell\bar{\nu}_{\ell})=\tau_{B_{s}}\times\int^{(m_{B_{s}}-m_{D_{s}})^{2}}_{0}\!\!\!dq^{2}\frac{d\Gamma(B_{s}\to D_{s}\ell\bar{\nu}_{\ell})}{dq^{2}},$$

(25)

where $q_{\rm max}^{2}=(m_{B_{s}}-m_{D_{s}})^{2}$ and $\tau_{B_{s}}$ is the $B_{s}$-meson lifetime. Here the differential decay widths is

	$\displaystyle\frac{d\Gamma(B_{s}\to D_{s}\ell\bar{\nu}_{\ell})}{dq^{2}}$	$\displaystyle=$	$\displaystyle\frac{G_{F}^{2}|V_{cb}|^{2}}{192\pi^{3}m_{B_{s}}^{3}}\left(1-\frac{m_{\ell}^{2}}{q^{2}}\right)^{2}\left[\left(1+\frac{m_{\ell}^{2}}{2q^{2}}\right)\right.$		(26)
			$\displaystyle\times\lambda^{\frac{3}{2}}(q^{2})|f_{+}^{B_{s}\to D_{s}}(q^{2})|^{2}$
			$\displaystyle+\frac{3m^{2}_{\ell}}{2q^{2}}(m^{2}_{B_{s}}-m^{2}_{D_{s}})^{2}$
			$\displaystyle\times\left.\lambda^{\frac{1}{2}}(q^{2})|f_{0}^{B_{s}\to D_{s}}(q^{2})|^{2}\right],$

where $\lambda(q^{2})=(m_{B_{s}}^{2}+m_{D_{s}}^{2}-q^{2})^{2}-4m_{B_{s}}^{2}m_{D_{s}}^{2}$, which is the phase-space factor. $|V_{cb}|$, $G_{F}$ and $m_{\ell}$ are CKM matrix element, the Fermi-coupling constant and the lepton mass, respectively, and we take Zyla:2020zbs : $\tau_{B_{s}}=(1.510\pm{0.004})\times 10^{-12}s$, $|V_{cb}|=(40.5\pm 1.5)\times 10^{-3}$, $G_{F}=1.1663787(6)\times 10^{-5}{\rm GeV}^{-2}$ and $\tau$-lepton mass $m_{\tau}=1.776\pm 0.00012{\rm GeV}$. For lepton $\ell^{\prime}=e$ or $\mu$, its mass is negligible, and then the above differential decay width can be simplified as

$$\frac{d\Gamma(B_{s}\to D_{s}\ell^{\prime}\bar{\nu}_{\ell^{\prime}})}{dq^{2}}=\frac{G_{F}^{2}|V_{cb}|^{2}}{192\pi^{3}m_{B_{s}}^{3}}\,\lambda^{3/2}(q^{2})|f_{+}^{B_{s}\to D_{s}}(q^{2})|^{2},$$

(27)

where $f_{0}^{B_{s}\to D_{s}}(q^{2})$ has zero contribution due to chiral suppression.

We present the differential decay widths of $B_{s}\to D_{s}\tau\bar{\nu}_{\tau}$ and $B_{s}\to D_{s}\ell^{\prime}\bar{\nu}_{\ell^{\prime}}$ in Figure 4, in which the solid lines are for the central choices of input parameters, and the shaded bands are uncertainties by adding all the errors caused by the error sources such as $f^{B_{s}\to D_{s}}_{+,0}(q^{2})$, $m_{B_{s}}$, $m_{D_{s}}$, $|V_{cb}|$, $G_{F}$ and $m_{\ell}$, etc., in quadrature. In addition, the predictions under the RQM approach Faustov:2012mt and the LQCD approach McLean:2019qcx are also given. One may observe that our prediction of $d\Gamma(B_{s}\to D_{s}\tau\bar{\nu}_{\tau})/dq^{2}$ is consistent with the LQCD and RQM predictions in Refs. McLean:2019qcx ; Faustov:2012mt ; And for $d\Gamma(B_{s}\to D_{s}\ell^{\prime}\bar{\nu}_{\ell^{\prime}})/dq^{2}$, our prediction agrees with the LQCD and RQM predictions McLean:2019qcx ; Faustov:2012mt in larger $q^{2}$ region, but is smaller than those predictions in lower $q^{2}$ region.