Studies of the resonance components in the $B_s$ decays into charmonia plus kaon pair

TL;DR

This paper investigates the resonance structures in $B_s$ meson decays to charmonia and kaon pairs using perturbative QCD, successfully modeling invariant mass spectra and predicting decay rates and polarization fractions relevant for future experiments.

Contribution

It provides a detailed analysis of resonance contributions in $B_s$ decays, including S, P, D waves, and predicts decay branching ratios and polarization fractions with comparison to LHCb data.

Findings

Invariant mass spectra well reproduced for most resonances.

Predicted decay branching ratios around 10^{-4}.

Polarization fractions align with existing LHCb measurements.

Abstract

In this work, the decays of $B_s$ mesons to a charmonium state and a $K^+K^-$ pair are carefully investigated in the perturbative QCD approach. Following the latest fit from the LHCb experiment, we restrict ourselves to the case where the produced $K^+K^-$ pair interact in isospin zero $S$, $P$, and $D$ wave resonances in the kinematically allowed mass window. Besides the dominant contributions of the $\phi(1020)$ resonance in the $P$-wave and $f_2'(1525)$ in the $D$-wave, other resonant structures in the high mass region as well as the $S$-wave components are also included. The invariant mass spectra for most of the resonances in the $B_s\rightarrow J/\psi K^+K^-$ decay are well reproduced. The obtained three-body decay branching ratios can reach the order of $10^{-4}$, which seem to be accessible in the near future experiments. The associated polarization fractions of those vector-vector and vector-tensor modes are also predicted, which are compared with the existing data from LHCb Collaboration.

Tables

Table 1. Table 1: The relevant resonance parameters in the B s → X ​ K + ​ K − → subscript 𝐵 𝑠 𝑋 superscript 𝐾 superscript 𝐾 B_{s}\rightarrow XK^{+}K^{-} decays.

Resonance	$J^{PC}$	Resonance formalism	Mass (MeV)	Width (MeV)	Source
$f_0\left(980\right)$	$0^{++}$	Flatté	990	$\cdots$	PDG pdg2018
$f_0\left(1370\right)$	$0^{++}$	BW	1475	113	LHCb prd86052006
$\varphi \left(1020\right)$	$1^{--}$	BW	1019	4.25	PDG pdg2018
$\varphi \left(1680\right)$	$1^{--}$	BW	1689	211	Belle prd80031101
$f_2\left(1270\right)$	$2^{++}$	BW	1276	187	PDG pdg2018
$f_2^{\prime }\left(1525\right)$	$2^{++}$	BW	1525	73	PDG pdg2018
$f_2\left(1750\right)$	$2^{++}$	BW	1737	151	Belle epjc32323
$f_2\left(1950\right)$	$2^{++}$	BW	1980	297	Belle epjc32323

Table 2. Table 2: Branching ratios of S 𝑆 S -wave resonance contributions to the B s → ( J / ψ , ψ ​ ( 2 ​ S ) , η c , η c ​ ( 2 ​ S ) ) ​ K + ​ K − → subscript 𝐵 𝑠 𝐽 𝜓 𝜓 2 𝑆 subscript 𝜂 𝑐 subscript 𝜂 𝑐 2 𝑆 superscript 𝐾 superscript 𝐾 B_{s}\rightarrow(J/\psi,\psi(2S),\eta_{c},\eta_{c}(2S))K^{+}K^{-} decays. Theoretical errors correspond to the uncertainties of Gegenbauer moments and hard scales, respectively.

Modes	$ℬ\left(R=f_0\left(980\right)\right)$	$ℬ\left(R=f_0\left(1370\right)\right)$	S-wave111We are not including contributions from the nonresonant $S$-wave.
$J/\psi K^{+}{K}^{-}$	$\left({4.3}_{-1.3-0.3}^{+1.3+0.7}\right)\times {10}^{-5}$	$\left({9.4}_{-3.4-0.4}^{+3.1+1.5}\right)\times {10}^{-6}$	$\left({3.4}_{-0.9-0.4}^{+1.0+0.5}\right)\times {10}^{-5}$
Data prd87072004	$\left(3.7\sim 9.2\right)\times {10}^{-5}$ 222We quote the range of measurement since the fit fraction of $f_0\left(980\right)$ is strongly parametrization dependent.	$\left({9.2}_{-9.2}^{+3.3}\right)\times {10}^{-6}$ 333The fit fraction statistical and systematic are added in quadrature.	$\cdots$
${\eta }_c{K}^{+}{K}^{-}$	$\left({4.6}_{-1.1-0.5}^{+0.8+0.8}\right)\times {10}^{-5}$	$\left({1.1}_{-0.3-0.1}^{+0.2+0.2}\right)\times {10}^{-5}$	$\left({3.5}_{-0.7-0.3}^{+0.8+0.6}\right)\times {10}^{-5}$
$\psi \left(2S\right)K^{+}{K}^{-}$	$\left({8.0}_{-2.6-1.2}^{+2.2+1.8}\right)\times {10}^{-6}$	$\left({8.3}_{-3.0-0.8}^{+3.0+1.9}\right)\times {10}^{-7}$	$\left({6.2}_{-1.7-0.8}^{+1.7+1.4}\right)\times {10}^{-6}$
${\eta }_c\left(2S\right)K^{+}{K}^{-}$	$\left({9.3}_{-1.9-1.1}^{+1.7+1.9}\right)\times {10}^{-6}$	$\left({1.4}_{-0.3-0.1}^{+0.3+0.3}\right)\times {10}^{-6}$	$\left({7.2}_{-1.3-0.9}^{+1.5+1.4}\right)\times {10}^{-6}$

Table 3. Table 3: Branching ratios of P 𝑃 P -wave resonance contributions to the B s → ( J / ψ , ψ ​ ( 2 ​ S ) , η c , η c ​ ( 2 ​ S ) ) ​ K + ​ K − → subscript 𝐵 𝑠 𝐽 𝜓 𝜓 2 𝑆 subscript 𝜂 𝑐 subscript 𝜂 𝑐 2 𝑆 superscript 𝐾 superscript 𝐾 B_{s}\rightarrow(J/\psi,\psi(2S),\eta_{c},\eta_{c}(2S))K^{+}K^{-} decays. Theoretical errors are attributed to the Gegenbauer moments, form factor ratios, and hard scales, respectively. The statistical and systematic uncertainties from data pdg2018 ; jhep070212017 ; jhep080372017 are combined in quadrature.

Modes	$ℬ\left(R=\varphi \left(1020\right)\right)$	$ℬ\left(R=\varphi \left(1680\right)\right)$	P-wave
$J/\psi K^{+}{K}^{-}$	${5.7}_{-0.2-0.7-0.0}^{+0.2+0.9+0.2}\times {10}^{-4}$	${3.6}_{-0.1-0.3-0.3}^{+0.1+0.4+0.3}\times {10}^{-5}$	${5.9}_{-0.1-0.7-0.0}^{+0.1+0.9+0.3}\times {10}^{-4}$
Data jhep080372017	$\left(5.6\pm 0.5\right)\times {10}^{-4}$ 111The fit fractions determined from the Dalitz plot analysis have been converted into the branching ratio measurements.	$\left(3.2\pm 0.4\right)\times {10}^{-5}$ 111The fit fractions determined from the Dalitz plot analysis have been converted into the branching ratio measurements.	$\cdots$
${\eta }_c{K}^{+}{K}^{-}$	${2.4}_{-0.2-0.4-0.0}^{+0.2+0.5+0.2}\times {10}^{-4}$	${1.1}_{-0.1-0.2-0.0}^{+0.1+0.2+0.0}\times {10}^{-5}$	${2.4}_{-0.1-0.3-0.0}^{+0.2+0.6+0.4}\times {10}^{-4}$
Data jhep070212017	$\left(2.5\pm 0.4\right)\times {10}^{-4}$ 222The experimental data is obtained by the product of $ℬ\left(B_s\to {\eta }_c\varphi \left(1020\right)\right)$ and $ℬ(\varphi \left(1020\right)\to K^{+}{K}^{-}$).	$\cdots$	$\cdots$
$\psi \left(2S\right)K^{+}{K}^{-}$	${2.4}_{-0.1-0.4-0.0}^{+0.0+0.3+0.0}\times {10}^{-4}$	${3.2}_{-0.1-0.4-0.0}^{+0.1+0.7+0.4}\times {10}^{-6}$	${2.3}_{-0.1-0.4-0.0}^{+0.1+0.3+0.1}\times {10}^{-4}$
Data pdg2018	$\left(2.6\pm 0.3\right)\times {10}^{-4}$ 333The experimental data is obtained by the product of $ℬ\left(B_s\to \psi \left(2S\right)\varphi \left(1020\right)\right)$ and $ℬ(\varphi \left(1020\right)\to K^{+}{K}^{-}$).	$\cdots$	$\cdots$
${\eta }_c\left(2S\right)K^{+}{K}^{-}$	${8.0}_{-0.6-1.8-0.5}^{+0.6+2.2+0.9}\times {10}^{-5}$	${1.1}_{-0.1-0.2-0.0}^{+0.1+0.2+0.2}\times {10}^{-6}$	${8.0}_{-0.7-1.8-0.5}^{+0.6+2.0+1.0}\times {10}^{-5}$

Table 4. Table 4: Branching ratios of D 𝐷 D -wave resonant contributions to the B s → ( J / ψ , ψ ​ ( 2 ​ S ) , η c , η c ​ ( 2 ​ S ) ) ​ K + ​ K − → subscript 𝐵 𝑠 𝐽 𝜓 𝜓 2 𝑆 subscript 𝜂 𝑐 subscript 𝜂 𝑐 2 𝑆 superscript 𝐾 superscript 𝐾 B_{s}\rightarrow(J/\psi,\psi(2S),\eta_{c},\eta_{c}(2S))K^{+}K^{-} decays. For theoretical errors, see Table. 3 . The statistical and systematic uncertainties from jhep080372017 are combined in quadrature.

Modes	$ℬ\left(R=f_2^{\prime }\left(1525\right)\right)$	$ℬ\left(R=f_2\left(1270\right)\right)$	$ℬ\left(R=f_2\left(1750\right)\right)$	$ℬ\left(R=f_2\left(1950\right)\right)$	D-wave
$J/\psi K^{+}{K}^{-}$	${8.9}_{-2.3-1.3-0.1}^{+2.8+1.8+0.9}\times {10}^{-5}$	${2.3}_{-0.5-0.4-0.1}^{+0.8+0.5+0.3}\times {10}^{-7}$	${5.0}_{-1.3-0.4-0.0}^{+1.5+0.6+0.1}\times {10}^{-6}$	${4.1}_{-1.0-0.7-0.0}^{+1.3+0.8+0.1}\times {10}^{-6}$	${9.2}_{-2.5-1.2-0.0}^{+2.8+1.9+0.9}\times {10}^{-5}$
Data jhep080372017	$\left(8.5\pm 1.2\right)\times {10}^{-5}$	$\left(1.3\pm 0.3\right)\times {10}^{-5}$	${4.7}_{-2.1}^{+2.4}\times {10}^{-6}$	${3.5}_{-1.4}^{+1.7}\times {10}^{-6}$	$\cdots$
${\eta }_c{K}^{+}{K}^{-}$	${4.9}_{-1.8-1.3-0.2}^{+2.1+1.6+0.5}\times {10}^{-5}$	${1.5}_{-0.5-0.3-0.1}^{+0.6+0.3+0.1}\times {10}^{-7}$	${2.3}_{-0.8-0.5-0.1}^{+1.0+0.8+0.3}\times {10}^{-6}$	${2.6}_{-0.9-0.4-0.0}^{+1.1+0.4+0.1}\times {10}^{-6}$	${4.9}_{-1.7-1.2-0.1}^{+2.2+1.6+0.4}\times {10}^{-5}$
$\psi \left(2S\right)K^{+}{K}^{-}$	${1.3}_{-0.4-0.2-0.1}^{+0.4+0.2+0.0}\times {10}^{-5}$	${5.9}_{-1.6-1.2-0.3}^{+1.8+1.4+0.3}\times {10}^{-8}$	$\cdots$	$\cdots$	${1.3}_{-0.4-0.2-0.1}^{+0.3+0.2+0.0}\times {10}^{-5}$
${\eta }_c\left(2S\right)K^{+}{K}^{-}$	${0.7}_{-0.2-0.2-0.1}^{+0.4+0.3+0.1}\times {10}^{-5}$	${3.2}_{-1.1-0.7-0.1}^{+1.4+0.8+0.1}\times {10}^{-8}$	$\cdots$	$\cdots$	${0.7}_{-0.2-0.2-0.0}^{+0.3+0.3+0.1}\times {10}^{-5}$

Table 5. Table 5: Polarization fractions for the decays B s → ( J / ψ , ψ ​ ( 2 ​ S ) ) ​ ϕ / f 2 ( → K + ​ K − ) → subscript 𝐵 𝑠 annotated 𝐽 𝜓 𝜓 2 𝑆 italic-ϕ subscript 𝑓 2 → absent superscript 𝐾 superscript 𝐾 B_{s}\rightarrow(J/\psi,\psi(2S))\phi/f_{2}(\rightarrow K^{+}K^{-}) . For theoretical errors, see Table.  3 . The experimental data are taken from Ref. jhep080372017 , where the statistical and systematic uncertainties are combined in quadrature.

Modes	$f_0(\%)$	$f_{\parallel }(\%)$	$f_{⟂}(\%)$
$B_s\to J/\psi \varphi \left(1020\right)\left(\to K^{+}{K}^{-}\right)$	${50.6}_{-1.5-4.4-1.7}^{+1.2+5.9+2.5}$	${24.4}_{-0.6-2.1-1.3}^{+1.1+1.8+1.0}$	${24.9}_{-0.5-3.6-1.1}^{+0.6+2.7+0.9}$
Data jhep080372017	$50.9\pm 0.4$	$23.1\pm 0.5$	$26.0\pm 0.6$
$B_s\to J/\psi \varphi \left(1680\right)\left(\to K^{+}{K}^{-}\right)$	${49.1}_{-1.1-3.6-0.9}^{+1.5+3.8+0.6}$	${20.1}_{-0.5-0.2-0.3}^{+0.2+0.3+0.7}$	${30.8}_{-1.1-3.6-0.3}^{+0.9+3.3+0.1}$
Data jhep080372017	$44.0\pm 3.9$	$32.7\pm 3.6$	$23.3\pm 3.6$
$B_s\to \psi \left(2S\right)\varphi \left(1020\right)\left(\to K^{+}{K}^{-}\right)$	${44.1}_{-2.4-5.0-1.8}^{+1.3+3.8+0.6}$	${23.2}_{-0.8-1.3-0.4}^{+1.2+1.8+1.1}$	${32.7}_{-0.5-2.4-0.2}^{+1.1+3.1+0.7}$
Data plb762253	$42.2\pm 1.4$	$\cdots$	$26.4\pm 2.4$
$B_s\to \psi \left(2S\right)\varphi \left(1680\right)\left(\to K^{+}{K}^{-}\right)$	${45.3}_{-0.4-3.3-1.3}^{+0.3+3.5+0.2}$	${15.1}_{-0.1-0.2-0.3}^{+0.2+0.3+0.9}$	${39.6}_{-0.1-3.8-0.0}^{+0.1+3.5+0.4}$
$B_s\to J/\psi f_2^{\prime }\left(1525\right)\left(\to K^{+}{K}^{-}\right)$	${51.1}_{-14.2-3.6-0.3}^{+14.1+4.3+1.6}$	${26.2}_{-7.6-2.1-0.4}^{+7.6+1.9+0.0}$	${22.7}_{-6.5-2.2-1.4}^{+6.7+1.7+0.6}$
Data jhep080372017	$46.8\pm 1.9$	$33.8\pm 2.3$	$19.4\pm 2.8$
$B_s\to J/\psi f_2\left(1270\right)\left(\to K^{+}{K}^{-}\right)$	${42.9}_{-14.2-2.7-2.4}^{+12.5+3.0+2.3}$	${29.5}_{-6.4-1.5-0.6}^{+7.4+1.3+0.6}$	${27.6}_{-6.1-1.6-1.7}^{+6.7+1.4+1.8}$
Data jhep080372017	$76.9\pm 5.5$	$6.0\pm 4.2$	$17.1\pm 5.0$
$B_s\to J/\psi f_2\left(1750\right)\left(\to K^{+}{K}^{-}\right)$	${53.7}_{-14.4-3.3-1.7}^{+13.9+1.8+2.1}$	${25.3}_{-7.6-0.7-0.3}^{+7.8+1.6+0.8}$	${21.0}_{-6.3-1.0-1.9}^{+6.5+1.7+0.9}$
Data jhep080372017	$58.2\pm 13.9$	$31.7\pm 12.4$	${10.1}_{-6.1}^{+16.8}$
$B_s\to J/\psi f_2\left(1950\right)\left(\to K^{+}{K}^{-}\right)$	${30.2}_{-11.1-1.2-0.0}^{+13.1+3.1+1.8}$	${36.9}_{-6.9-1.5-1.4}^{+5.8+0.6+0.6}$	${32.9}_{-6.2-1.6-2.5}^{+5.2+0.6+0.8}$
Data jhep080372017	${2.2}_{-1.5}^{+6.7}$	$38.3\pm 13.8$	$59.5\pm 14.2$
$B_s\to \psi \left(2S\right)f_2^{\prime }\left(1525\right)\left(\to K^{+}{K}^{-}\right)$	${41.9}_{-13.5-2.7-1.9}^{+14.2+1.8+1.0}$	${34.7}_{-8.5-2.2-0.0}^{+8.0+2.8+0.6}$	${23.4}_{-5.7-0.2-1.5}^{+5.5+0.5+1.5}$
$B_s\to \psi \left(2S\right)f_2\left(1270\right)\left(\to K^{+}{K}^{-}\right)$	${36.4}_{-12.4-2.1-1.0}^{+13.9+3.2+0.9}$	${37.0}_{-8.1-1.7-0.6}^{+7.3+1.1+0.8}$	${26.5}_{-5.8-1.4-1.6}^{+5.2+1.0+1.7}$

Equations

$$B_s(p_B)\to X(p_3)K^{+}(p_1)K^{-}(p_2).$$

$$p_1=(\zeta p^{+},\eta (1-\zeta )p^{+},{\text{P}}_{\text{T}}),p_2=((1-\zeta )p^{+},\eta \zeta p^{+},-{\text{P}}_{\text{T}}),p_3=\frac{M}{\sqrt{2}}(r^2,1-\eta ,{\text{0}}_{\text{T}}),$$

$$k_B$$

$${\Phi }_S^{I=0}=\frac{1}{\sqrt{2N_c}}[to0.0pt/\text{hss}p{\varphi }_S^0(z,\zeta ,{\omega }^2)+\omega {\varphi }_S^s(z,\zeta ,{\omega }^2)+\omega (to0.0pt/\text{hss}nto0.0pt/\text{hss}v-1){\varphi }_S^t(z,\zeta ,{\omega }^2)],$$

$${\varphi }_S^0(z,\zeta ,{\omega }^2)$$

$${\varphi }_S^s(z,\zeta ,{\omega }^2)$$

$${\varphi }_S^t(z,\zeta ,{\omega }^2)$$

$$F_S({\omega }^2)=[\frac{m_{f_0(980)}^2}{m_{f_0(980)}^2-{\omega }^2-i{m}_{f_0(980)}(g_{\pi \pi }{\rho }_{\pi \pi }+g_{KK}{\rho }_{KK}{F}_{KK}^2)}+c_{f_0(1370)}B{W}_{f_0(1370)}({\omega }^2)](1+c_{f_0(1370)}{)}^{-1}.$$

$$B{W}_R({\omega }^2)=\frac{m_R^2}{m_R^2-{\omega }^2-i{m}_R\Gamma ({\omega }^2)},$$

$$\Gamma ({\omega }^2)={\Gamma }_0(\frac{q_K}{q_{K0}}{)}^{2L_R+1}\frac{m_R}{\omega }{F}_R^2.$$

$$F_R=\{\begin{array}{rlr}& 1& L_R=0,\\ & \frac{\sqrt{1+z_0}}{\sqrt{1+z}}& L_R=1,\\ & \frac{\sqrt{z_0^2+3z_0+9}}{\sqrt{z^2+3z+9}}& L_R=2,\end{array}$$

$${\Phi }_P^L$$

$${\Phi }_P^T$$

$${\varphi }_P^0(z,\zeta ,\omega )$$

$${\varphi }_P^s(z,\zeta ,\omega )$$

$${\varphi }_P^t(z,\zeta ,\omega )$$

$${\varphi }_P^T(z,\zeta ,\omega )$$

$${\varphi }_P^a(z,\zeta ,\omega )$$

$${\varphi }_P^v(z,\zeta ,\omega )$$

$$F_P^{\parallel }({\omega }^2)=[B{W}_{\varphi (1020)}({\omega }^2)+c_{\varphi (1680)}B{W}_{\varphi (1680)}({\omega }^2)](1+c_{\varphi (1680)}{)}^{-1}.$$

$${\Phi }_D^L$$

$${\Phi }_D^T$$

$${\varphi }_D^0(z,\zeta ,\omega )$$

$${\varphi }_D^T(z,\zeta ,\omega )$$

$$F_D^{\parallel }({\omega }^2)=\sum _i{c}_{i}B{W}_i({\omega }^2),$$

$${\varphi }_D^s(z,\zeta ,\omega )$$

$${\varphi }_D^t(z,\zeta ,\omega )$$

$${\varphi }_D^a(z,\zeta ,\omega )$$

$${\varphi }_D^v(z,\zeta ,\omega )$$

$$⟨f_2(p,\lambda )|j_{\mu \nu }(0)|0⟩$$

$$⟨f_2(p,\lambda )|j_{\mu \nu \rho }(0)|0⟩$$

$$⟨K^{+}(p_1)K^{-}(p_2)|j_{\mu \nu }(0)/j_{\mu \nu \rho }(0)|0⟩\approx \sum _{\lambda }⟨K^{+}(p_1)K^{-}(p_2)|f_2(p,\lambda )⟩\frac{1}{{\mathcal{D}}_{f_2}}⟨f_2(p,\lambda )|j_{\mu \nu }/j_{\mu \nu \rho }|0⟩,$$

$$⟨K^{+}(p_1)K^{-}(p_2)|j_{\mu \nu }|0⟩$$

$$⟨K^{+}(p_1)K^{-}(p_2)|j_{\mu \nu \rho }|0⟩$$

$$q_{\mu }{q}_{\nu }-\frac{1}{3}{p}_{\mu }{p}_{\nu }=\frac{2}{3}(1-6\zeta +6{\zeta }^2)p_{\mu }{p}_{\nu },$$

$$(q_{\mu }{p}_{\nu }-q_{\nu }{p}_{\mu })q_{\rho }=2(2\zeta -1)\sqrt{\zeta (1-\zeta )}{\omega }^3,$$

$$\mathcal{L}(\zeta )=1-6\zeta +6{\zeta }^2,\mathcal{T}(\zeta )=(2\zeta -1)\sqrt{\zeta (1-\zeta )},$$

$$d_{00}^{l=2}(\theta )$$

Full text

I Introduction

The three-body mode $B_{s}\rightarrow J/\psi K^{+}K^{-}$ is of particular interest in searches for intermediate states in the $B_{s}$ decay chain. Since the LHCb Collaboration prd87072004 found no obvious structures in the $J/\psi K^{+}$ invariant mass distribution, the $B_{s}\rightarrow J/\psi K^{+}K^{-}$ decay proceeds predominantly via $B_{s}\rightarrow J/\psi R$ with the quasi-two-body intermediate state $R$ subsequently decaying into $K^{+}K^{-}$. For the concerned $B_{s}$ decay, the $K^{+}K^{-}$ system arise from pure $s\bar{s}$ source, and thus these resonances are isoscalar. Taking into account the conservation of $P$-parity and $C$-parity, the produced resonances are limited to quantum numbers $J^{PC}=0^{++},1^{--},2^{++},...$ with isospin $I=0$. Among them, the largest component comes from the $\phi(1020)$ in a $P$-wave configuration pdg2018 . Several Collaborations prd87072004 ; prd546596 ; prd88114006 have presented a measurement of the $B_{s}\rightarrow J/\psi\phi(1020)$ mode with $\phi(1020)$ decays to $K^{+}K^{-}$. The current world averages of the absolute branching ratio $\mathcal{B}(B_{s}\rightarrow J/\psi\phi(1020))$ can reach the order of $10^{-3}$. Another $P$-wave resonance $\phi(1680)$, whose contribution has more than 2 statistical standard deviation $(\sigma)$ significance, is also included by the LHCb Collaboration prd87072004 in its best fit model. Two well known scalar resonances, the $f_{0}(980)$ and $f_{0}(1370)$, are observed in the $K^{+}K^{-}$ mass spectrum by LHCb prd87072004 , which is the only data set available so far for the $S$-wave resonant structures.

Contributions from $D$-wave resonances are known to be non-negligible in this decay. The first observation of the decay sequence $B_{s}\rightarrow J/\psi f_{2}^{\prime}(1525),f_{2}^{\prime}(1525)\rightarrow K^{+}K^{-}$, was recently reported by the LHCb Collaboration prl108151801 , and later confirmed by the D0 Collaboration prd86092011 . Subsequently, the LHCb Collaboration prd87072004 have determined the final state composition of the decay channel using a modified Dalitz plot analysis where the decay angular distributions are included. The best fit model includes a nonresonant component and eight resonance states, whose absolute branching ratios are measured relative to that of the normalization decay mode $B^{+}\rightarrow J/\psi K^{+}$. In contrast to hadron collider experiments, the Belle Collaboration prd88114006 normalize to the absolute number of $B^{0}_{s}\bar{B}^{0}_{s}$ pairs produced and also present a measurement of the entire $B_{s}\rightarrow J/\psi K^{+}K^{-}$ components including resonant and nonresonant decays. More recently, the LHCb Collaboration jhep080372017 improved their measurements, in which the fit fractions of six resonances including $\phi(1020),\phi(1680),f_{2}(1270),f_{2}^{\prime}(1525),f_{2}(1750),f_{2}(1950)$ together with a $S$-wave structure in $B_{s}\rightarrow J/\psi K^{+}K^{-}$ are determined.

Above measurements have caught theoretical attention recently. The three-body decay $B_{s}\rightarrow J/\psi K^{+}K^{-}$ including its dominant contributions of the resonances $\phi(1020)$ and $f_{2}^{\prime}(1525)$ have been studied prd89095026 ; prd95036013 and the associated branching ratios have been obtained based on the framework of the factorization approach. Some recent analyses prd90094006 ; plb73770 ; prd90114004 had been carried out for $B/B_{s}$ decays into $J/\psi$ and the scalar, vector, and tensor resonances using chiral unitary theory, for which these states are shown to be generated from the meson-meson interaction. In Refs prd79074024 ; jhep09074 ; jhep02009 , the $K^{+}K^{-}$ $S$-wave contribution in the $\phi(1020)$ resonance region is estimated to be of the order $1-10\%$, in agreement with previous measurements from LHCb prd87072004 ; LHCb2012002 , CDF prl109171802 , and ATLAS jhep12072 . The significant $S$-wave effects may affect measurements of the $CP$ violating phase $\beta_{s}$ prd79074024 ; jhep09074 ; prd82076006 .

In this paper, we will consider the three-body $B_{s}$ decays involving charmonia and kaon pair in the final state under the quasi-two-body approximation in the framework of perturbative QCD approach (PQCD) prl744388 ; plb348597 . The factorization formalism for the three-body decays can be simplified to that for the two-body cases with the introduction of two-kaon distribution amplitudes (DAs), which absorb the strong interaction related to the production of the two kaon system. For the detailed description of the three-body nonleptonic $B$ decays in this approach, one can refer to plb561258 ; prd70054006 . The PQCD approach so far, has been successfully applied to the studies of the resonance contributions to the three-body $B/B_{s}$ decays in several recent papers prd91094024 ; plb76329 ; plb788468 ; epjc77199 ; prd97033006 ; prd98113003 ; 181112738 ; prd95056008 ; prd96036014 ; prd98056019 ; prd96093011 . As advanced before, the decays under study are dominated by a series of resonances in $S$, $P$ and $D$ waves, while contributions from resonances with spin greater than two are not expected since they are well beyond the available phase space. Each partial wave contribution is parametrized into the corresponding timelike form factors involved in the two-kaon DAs. For each partial wave form factor, we adopt the form as a linear combination of those resonances with the same spin. In the present paper, we take into account the following resonances 111 In the following, we also use the abbreviation $f_{0}$, $\phi$, and $f_{2}$ to denote the $S$, $P$, and $D$-wave resonances for simplicity.: $f_{0}(980)$, $f_{0}(1370)$, $\phi(1020)$, $\phi(1680)$, $f_{2}(1270)$, $f^{\prime}_{2}(1525)$, $f_{2}(1750)$, $f_{2}(1950)$, two scalar, two vector, and four tensor resonances in the context of the data presented in Refs prd87072004 ; jhep080372017 . All resonances are commonly described by Breit Wigner (BW) distributions, except for the $f_{0}(980)$ state, which is modelled by a Flatté function plb63228 .

Our presentation is divided as follows. In Sec. II, we present our model kinematics and describe the two-kaon DAs in different partial waves. The calculated branching ratios and polarizations for each resonance in the considered three-body decays as well as the numerical discussions are presented in Sec. III, and finally, conclusions are drawn in Sec IV. The factorization formulas for the decay amplitudes are collected in the Appendix.

II Kinematics and the two-kaon distribution amplitudes

Let us begin with the definition of the kinematic variables. It is convenient to work in the light-cone coordinates for the four-momenta of the initial and final states. The momentum of the decaying $B_{s}$ meson in its rest frame is chosen as $p_{B}=\frac{M}{\sqrt{2}}(1,1,\textbf{0}_{T})$ with the $B_{s}$ meson mass $M$. The momenta of the decay products will be denoted as $p_{1},p_{2}$ for the two kaons, and $p_{3}$ for the charmonia, with the specific charge assignment according to

[TABLE]

The momenta of three final states are defined as

[TABLE]

where $\zeta=p_{1}^{+}/p^{+}$ with $p^{+}=M(1-r^{2})/\sqrt{2}$ is the kaon momentum fraction. The mass ratio $r=m/M$ with the charmonium mass $m$. The kaon transverse momenta is expressed as $\textbf{P}_{\text{T}}=(\omega\sqrt{\zeta(1-\zeta)},0)$. The factor $\eta=p^{2}/(M^{2}-m^{2})$ is defined in terms of the invariant mass squared of the kaon pair $p^{2}=\omega^{2}$, which satisfies the momentum conservation $p=p_{1}+p_{2}=p_{B}-p_{3}$. The valence quark momenta labeled by $k_{B}$, $k_{3}$, and $k$, as indicated in Fig. 1 (a), are parametrized as

[TABLE]

in which $x_{B}$, $x_{3}$, $z$ denote the longitudinal momentum fractions, and $k_{iT}$ represent the transverse momenta. Since the light spectator quark momentum $k$ moves with the kaon pair in the plus direction, the minus component of its parton momentum should be very small, thus it can be neglected in the hard kernel, and then integrated out in the definition of the two-kaon distribution amplitudes. We also dropped $k_{B}^{+}$ because it vanishes in the hard amplitudes.

Since the $B_{s}$ meson wave function and the charmonium distribution amplitudes have successfully described various hadronic two-body and three-body charmonium $B$ decays epjc77199 ; prd97033006 ; prd98113003 ; 181112738 ; epjc77610 , we use the same ansatz as them. For the sake of brevity, their explicit expressions are not shown here and can be found in Refs. epjc77610 ; prd90114030 ; epjc75293 . Below, we briefly describe the two-kaon DAs in three partial waves and the associated form factors.

II.1 $S$-wave two-kaon DAs

The $S$-wave two-kaon DAs are introduced in analogy with the case of two-pion ones prd91094024 ; plb561258 , which are organized into

[TABLE]

with the null vectors $n=(1,0,\textbf{0}_{\text{T}})$ and $v=(0,1,\textbf{0}_{\text{T}})$. In what follows the subscripts $S$, $P$, and $D$ always associate with the corresponding partial waves. Above various twists DAs have similar forms as the corresponding twists for a scalar meson by replacing the scalar decay constant with the scalar form factor plb730336 , we adopt their asymptotic models as shown below prd91094024 ; 181112738 :

[TABLE]

with the isoscalar scalar form factor $F_{S}(\omega^{2})$ and the Gegenbauer moment $a_{1}$. Bearing in mind that only odd moments contribute in case of neutral scalar resonances owing to charge conjugation invariance or conservation of vector current prd73014017 . Therefore the first term in leading twist DAs come from $a_{1}$. Since the coefficients in the Gegenbauer expansion of the dimeson DAs are poorly known, we limit ourselves to leading term in the expansion.

For the scalar resonances, we include here only the components $f_{0}(980)$ and $f_{0}(1370)$, which are well established in the best fit model by the LHCb Collaboration prd87072004 . For the former, we use a Flatté description, while the latter is modelled by BW functions. The scalar form factor $F_{S}(\omega^{2})$ can be written as

[TABLE]

Hereafter, $c_{R}$ refers to the weight coefficient of the resonance $R$, to be determined by data. Their values are given in the next section. In what follows, all resonances $R$ with different quantum numbers will be labeled by the single letter $R$, without pointing to its quantum numbers. The two phase-space factors are $\rho_{\pi\pi}=2q_{\pi}/\omega$, $\rho_{KK}=2q_{K}/\omega$, where $q_{\pi(K)}$ is the pion (kaon) momentum in the dipion (dikaon) rest frame. The exponential factor $F_{KK}=e^{-\alpha q_{K}^{2}}$ with $\alpha=2.0\pm 0.25$ GeV*-2* prd78074023 ; prd89092006 is introduced above the $KK$ threshold and serves to reduce the $\rho_{KK}$ factor as the invariant mass increases. The constants $g_{\pi\pi}$ and $g_{KK}$ are the $f_{0}(980)$ couplings to $\pi\pi$ and $K\bar{K}$ final states respectively. We use $g_{\pi\pi}=167$ MeV and $g_{KK}/g_{\pi\pi}=3.47$ as determined by LHCb prd90012003 . The BW amplitude in generic form is

[TABLE]

where $m_{R}$ is the resonance pole mass and $\Gamma(\omega^{2})$ is its energy-dependent width which may be parametrized in a form that ensures the correct behavior near threshold,

[TABLE]

Here $\Gamma_{0}$ and $q_{K0}$ are $\Gamma(\omega^{2})$ and $q_{K}$, evaluated at the resonance pole mass, respectively. $L_{R}$ is the orbital angular momentum in the $K^{+}K^{-}$ decay and is equal to the spin of resonance $R$ because kaons have spin 0. The $L_{R}=0,1,2,...$ correspond to the $S,P,D,...$ partial wave resonances. The Blatt-Weisskopf barrier factors $F_{R}$ TNP for scalar, vector and tensor states are

[TABLE]

with $z=r^{2}q_{K}^{2}$ and $z_{0}$ represents the value of $z$ when $\omega=m_{R}$. The meson radius parameters $r$ are dependent on the momentum of the decay particles in the parent rest frame. Modifying the $r$ changes slightly our results, as discussed in the next section. Hence, we set this parameter to be 1.5 $\text{GeV}^{-1}$ (corresponding to $0.3$ fm) for all the considered resonances, as is obtained in prd87072004 .

II.2 $P$-wave two-kaon DAs

In Ref. prd98113003 we have constructed the $P$-wave DAs including both longitudinal and transverse polarizations for the pion pair. Naively, the $P$-wave two-kaon ones can be obtained by replacing the pion vector form factors by the corresponding kaon ones. The explicit expressions read

[TABLE]

where the superscripts $L$ and $T$ on the left-hand side denote the longitudinal polarization and transverse polarization, respectively. Here $\epsilon^{\mu\nu\rho\sigma}$ is the totally antisymmetric unit Levi-Civita tensor with the convention $\epsilon^{0123}=1$. The transverse polarization vector $\epsilon_{T}$ for the dikaon system has the same form as that of dipion prd98113003 . The various twists DAs in Eq. (II.2) can be expanded in terms of the Gegenbauer polynomials:

[TABLE]

where the two $P$-wave form factors $F_{P}^{\parallel}$ and $F_{P}^{\perp}$ serve as the normalization of the two-kaon DAs. They play a similar role with the vector and tensor decay constants in the definition of the vector meson DAs prd65014007 . The Gegenbauer moments $a_{2}^{i}$ will be regarded as free parameters and determined in this work.

As mentioned in the Introduction, the form factor $F^{\parallel}_{P}$ is given by the coherence summation of the two vector resonances $\phi(1020)$ and $\phi(1680)$,

[TABLE]

According to the argument in plb76329 [see Eq.(12)], motivated by the pole model, each form factor is proportional to the decay constant associated with each resonance state. Therefore, for the $F^{\perp}_{P}$ , we assume it have the same phase as $F^{\parallel}_{P}$ and employ the approximate relation $F_{P}^{\perp}/F_{P}^{\parallel}\sim f^{T}_{V}/f_{V}$ with $f^{T}_{V}$ ($f_{V}$) being the tensor (vector) decay constant for the corresponding vector meson in the following calculations. It is worth stressing that the current data are still not sufficient to determine the two form factors separately. In principle, the longitudinal decay constant $f_{V}$ could be extracted from the measurements, while the transverse one $f^{T}_{V}$ has to be calculated in the QCD sum rule or Lattice QCD technique. It should be noted that the latter one is renormalization scheme dependent and renormalization scale dependent, respectively. Then the approximate relation $F_{P}^{\perp}/F_{P}^{\parallel}$ could vary with the choices of decay constants at different energy scales. In the numerical analysis, their values are chosen at the typical scale $\mu=1$ GeV, which enters the perturbative calculation in PQCD.

II.3 $D$-wave two-kaon DAs

Recalling that the tensor meson DAs are constructed in analogy with the vector ones by introducing a new polarization vector $\epsilon_{\bullet}$ which is related to the polarization tensor $\epsilon_{\mu\nu}(\lambda)$ with helicity $\lambda$ prd83014008 . Following a similar procedure, we decompose the $D$-wave two-kaon DAs associated with longitudinal and transverse polarizations into

[TABLE]

respectively, where the prefactor $\sqrt{\frac{2}{3}}$ ($\sqrt{\frac{1}{2}}$) comes from the different definitions of the polarization vector between the vector and tensor mesons for the longitudinal (transverse) polarization. The leading twist DAs $\Phi_{D}^{0}$ and $\Phi_{D}^{T}$ have similar asymptotic forms as the corresponding ones for a tensor meson. More precisely,

[TABLE]

Here, we solely employ the first nonvanishing leading term in the expansion for previously mentioned reasons. The moments $a_{1}^{0}$ and $a_{1}^{T}$ are regarded as free parameters and determined in the next section. Note that the $\zeta$ dependent terms $\mathcal{L}(\zeta)$ and $\mathcal{T}(\zeta)$ are different from those in Eq. (II.2). We will derive their expressions later. The kaon tensor form factor $F^{\parallel}_{D}(\omega^{2})$ can be represented by

[TABLE]

where the summation is performed over the intermediate tensor mesons: $f^{\prime}_{2}(1525)$, $f_{2}(1270)$, $f_{2}(1750)$, $f_{2}(1950)$. $c_{i}$ are the corresponding weight coefficients. The expressions for the twist-3 DAs can be derived through the Wandzura-Wilczek relations as prd82054019 ; prd83034001

[TABLE]

Next we derive the $\zeta$ dependent terms for both the longitudinal and transverse polarization DAs. The two decay constants $f_{T}$ and $f_{T}^{T}$ of a tensor meson are defined by sandwiching the corresponding local current operators between the vacuum and a tensor meson prd82054019 ; prd83014008

[TABLE]

where $p$ and $m_{T}$ are momentum and mass of the tensor meson, respectively. The two interpolating currents $j_{\mu\nu}(0)$ and $j_{\mu\nu\rho}(0)$ are defined in prd82054019 ; prd83014008 . Let us begin with the local matrix element $\langle K^{+}(p_{1})K^{-}(p_{2})|j_{\mu\nu}(0)/j_{\mu\nu\rho}(0)|0\rangle$ associated with the $D$-wave form factors. Under the tensor-meson-dominant hypothesis prd82054019 , inserting the tensor intermediate in above matrix element, we get

[TABLE]

with $\mathcal{D}_{f_{2}}$ the resonance propagator 160503889 . The coupling constant $g_{f_{2}KK}$ is defined by the matrix element $\langle K^{+}(p_{1})K^{-}(p_{2})|f_{2}(p,\lambda)\rangle=\frac{g_{f_{2}KK}}{m_{T}}\epsilon_{\mu\nu}(\lambda)q^{\mu}q^{\nu}$ with $q=p_{1}-p_{2}$ prd82054019 . When applying the formula Eq. (II.3) and the completeness relation $\sum_{\lambda}\epsilon_{\mu\nu}(\lambda)\epsilon^{*}_{\rho\sigma}(\lambda)=\frac{1}{2}M_{\mu\rho}M_{\nu\sigma}+\frac{1}{2}M_{\mu\sigma}M_{\nu\rho}-\frac{1}{3}M_{\mu\nu}M_{\rho\sigma}$ with $M_{\mu\nu}=g_{\mu\nu}-p_{\mu}p_{\nu}/m^{2}_{T}$, Eq. (18) then leads to the following equations explicitly

[TABLE]

Note that the last terms above are power suppressed and can be omitted, because in our power counting, a light hadron mass is counted as a low scale relative to the heavy quark mass. Utilizing the approximation relation $q_{\mu}=(p_{1}-p_{2})_{\mu}\approx(2\zeta-1)p_{\mu}$ prd98113003 , one get for Eq. (19)

[TABLE]

in which the coefficient $1-6\zeta+6\zeta^{2}$ is absorbed into the longitudinal polarization DAs, giving rise to its $\zeta$ dependence. The matrix element in Eq. (20) for the choice $\mu,\nu,\rho=+,-,x$ is proportional to

[TABLE]

where the kinematic variables in Eq. (2) are used. Note that the contributions from other possible combination of the three Lorentz indexes are either zero or power suppressed. Then the $\zeta$ dependent factors of the longitudinal and transverse polarization DAs can be written as

[TABLE]

respectively. Above expressions can also be checked from the partial wave expansions of the helicity amplitudes. As is well known, the helicity 0 component is expanded in terms of Legendre polynomials $P_{l}(\text{cos}\theta)$, while the helicity $\pm 1$ ones proceeding in derivatives of the Legendre polynomials $P^{\prime}_{l}(\text{cos}\theta)$ jhep02009 . For $l=2$ $D$-wave amplitudes, the helicity angle $\theta$ is encoded into the Wigner-d functions, schematically pdg2018 :

[TABLE]

Following a similar prescription in the dipion system npb905373 ; prd96051901 , the polar angle $\theta$ of the $K^{+}$ in the rest frame of the dikaon is related to $\zeta$ through the relation

[TABLE]

with $m_{K}$ the kaon mass. Neglecting the kaon mass and employing Eq. (II.3), we also arrive at Eq. (23).

II.4 The differential branching ratio

The double differential branching ratio reads pdg2018

[TABLE]

where the differential variable $d\text{cos}\theta$ is replaced by $d\zeta$ via Eq. (25) in the limit of massless. The three-momenta of the kaon and charmonium in the rest reference frame of the $K\bar{K}$ system are given by

[TABLE]

respectively, with the standard K$\ddot{a}$ll$\acute{e}$n function $\lambda(a,b,c)=a^{2}+b^{2}+c^{2}-2(ab+ac+bc)$. The complete amplitude $\mathcal{A}$ through resonance intermediate for the concerned decay is decomposed into

[TABLE]

where $\mathcal{A}_{S}$, $\mathcal{A}_{P}$, and $\mathcal{A}_{D}$ denote the corresponding three partial wave decay amplitudes. Since the $\zeta$-dependent terms appear as an overall factor in each partial wave decay amplitudes, integrating the double differential distribution of Eq. (26) over $\zeta$ gives for the differential invariant mass branching ratio

[TABLE]

where the factors $1/3,1/5,\cdots$ are extracted from the individual helicity amplitudes for the integral of $\zeta$. The terms $\mathcal{A}^{0}$, $\mathcal{A}^{\parallel}$, and $\mathcal{A}^{\perp}$ represent the longitudinal, parallel, and perpendicular polarization amplitudes in the transversity basis, respectively. 222 The last two terms do not appear for those modes involving spinless $\eta_{c}/\eta_{c}(2S)$ in the final state. Note that interference between different partial wave vanishes because the $\zeta$ functions in Eqs. (II.1), (II.2), and (23), corresponding to $S$, $P$, and $D$ partial waves, are orthogonal. In the $J/\psi$ and $\psi(2S)$ cases, the decay amplitudes $\mathcal{A}_{S}$ and $\mathcal{A}_{P}$ here can be straightforwardly obtained from the previous publications prd91094024 ; prd98113003 by replacing the two-pion form factors and all pion masses and momenta with the respective kaon quantities. For $\mathcal{A}_{D}$, its factorization formula can be related to $\mathcal{A}_{P}$ by making the following replacement,

[TABLE]

For the involved $\eta_{c}$ and $\eta_{c}(2S)$ modes, the partial wave decay amplitudes are provided in Appendix.

III Numerical results

We first summarize all parameter values required for numerical applications. For the masses appearing in $B_{s}$ decays, we shall use the following values (in units of GeV) pdg2018 :

[TABLE]

The information on the decay constants (in units of GeV), the Wolfenstein parameters, together with the lifetime of $B_{s}$ mesons are adopted as prd82054019 ; prd98113003 ; epjc77610 ; prd90114030 ; epjc75293

[TABLE]

The masses and widths of the BW resonances are listed in Table 1, while the Flatté parameters for the $f_{0}(980)$ have been given in the previous section.

As mentioned before, the form factor ratio $r^{T}(R)=F^{\perp}/F^{\parallel}$ is approximately equal to the ratio of two decay constants $f^{T}/f$. From the numbers in Eq. (III), we have

[TABLE]

Whereas for other high states, since their decay constants are not known yet, we treat them as free parameters.

Now, we collect all the phenomenologically motivated parameters, such as weight coefficients $c_{i}$, Gegenbauer moments $a_{i}$, and form factor ratios $r^{T}(R)$, in each partial wave. Their central values are fixed to be

[TABLE]

When fitting to the experimental data, we assume that the concerned resonances with the same spin share the same set of the Gegenbauer moments. For the $S$-wave sector, the two experimental results of the $f_{0}(980)$ and $f_{0}(1370)$ components in Ref. prd87072004 are used to fit out $a_{1}$ and $c_{f_{0}(1370)}$. In Fig. 2 (a) and (b), we show the dependence of the branching ratios of the $f_{0}(980)$ and $f_{0}(1370)$ components as well as their combination in the $B_{s}\rightarrow J/\psi K^{+}K^{-}$ decay on the Gegenbauer moment $a_{1}$ and the phase of $c_{f_{0}(1370)}$, respectively. The module $|c_{f_{0}(1370)}|$ is chosen as 0.12 to maximize the overlap between the predicted curves and the experimental range. Apparently, both the $f_{0}(980)$ and $f_{0}(1370)$ modes can meet the data as setting $a_{1}\sim 0.8$ in Fig. 2 (a). This value is much larger than the corresponding $a_{1}=0.2$ that obtained in Ref. prd91094024 . The discrepancy is understandable with respect to the different nonperturbative dynamics of $f_{0}(980)$ decaying to the $KK$ and $\pi\pi$ pairs. In Fig. 2 (b), it is reflected that the constructive or destructive interference pattern between the two resonances vary with the phase. Here its value is taken to be $-\pi/2$ since the LHCb’s data prd87072004 favor the destructive interference. Of course, considering the sizeable experimental uncertainties, these parameters are difficult to be restricted precisely at this moment. As a case study with rough estimation, the related treatment about these parameters in this work is just a try. A convincible research should be performed through a global fit to more rich measurements in the future.

For the $P$-wave ones, we first use the experimental branching ratios of three decay channels $B_{s}\rightarrow J/\psi\phi(1020)(\rightarrow K^{+}K^{-})$ (longitudinal) jhep080372017 , $B_{s}\rightarrow\psi(2S)\phi(1020)(\rightarrow K^{+}K^{-})$ (longitudinal) pdg2018 , and $B_{s}\rightarrow\eta_{c}\phi(1020)(\rightarrow K^{+}K^{-})$ jhep070212017 to fit the three longitudinal Gegenbauer moments $a_{2}^{0}$, $a_{2}^{s}$, and $a_{2}^{t}$, then the three transverse ones can be constrained by the transverse polarization fractions of the former two modes. Finally, according to the fit fraction and polarizations of the $\phi(1680)$ component in the $B_{s}\rightarrow J/\psi K^{+}K^{-}$ decay from the LHCb jhep080372017 , one can determine the values of $r^{T}(\phi(1680))$ and $c_{\phi(1680)}$.

Since the $f^{\prime}_{2}(1525)$ component in the $J/\psi$ mode is well measured with a relatively high accuracy comparing with other $D$-wave resonances by the LHCb Collaboration jhep080372017 , we can exactly determine its weight coefficient $c_{f_{2}^{\prime}(1525)}$ and two Gegenbauer moments $a_{1}^{0}$ and $a_{1}^{T}$ based on its fit fraction and three polarizations. Following a similar procedure as above, we can determine $r^{T}(R)$ and the module of $c_{R}$ for other tensor resonances. As pointed out in epjc3941 , the form factor $F(s)$ with the time-like momentum transfer squared $s>4m_{K}^{2}$ could be analytically continued to the space-like region $s<0$. It has been known that a form factor is normalized to unity at $s=0$, because a soft probe cannot reveal the structure of a bound state. Therefore, we postulate that the kaon form factors should be constrained by such normalization condition. According to our fitted modules of the $D$-wave weight coefficients in Eq. (III), the phases of $c_{f_{2}(1270)}$ and $c_{f_{2}(1750)}$ are set to $\pi$ to ensure the normalization of the form factor $F^{\parallel}_{D}(0)=1$. Strictly speaking, the phases of the various coefficient $c_{i}$ in Eq. (III) should be determined from the interference fit fractions. However, the current available data are not yet sufficiently precise to extract them. Furthermore, the $f^{\prime}_{2}(1525)$ dominates over the $D$-wave contributions as shown below, the relative phases among these $c$ parameters have little effect on the total $D$-wave decay branching ratios, and our choice of these phases does not affect the magnitude estimation of either the individual resonance or the total contribution.

The calculated branching ratios of $S$, $P$, and $D$-wave resonance contributions to the $B_{s}\rightarrow(J/\psi,\psi(2S),\eta_{c},\eta_{c}(2S))K^{+}K^{-}$ decays are collected in Tables 2, 3, and 4, respectively. The last column of each Table are the corresponding total partial wave branching ratios. The theoretical errors stem from the uncertainties for fitted values of Gegenbauer moments $a_{i}$, the form factor ratios $r^{T}(R)$, and the hard scales $t$, respectively. For Gegenbauer moments in the twist-2 DAs, we vary their values within a $20\%$ range for the error estimation. The uncertainty of the ratio $r^{T}(R)$ in Eqs. (33) and (III) are general assigned to be $\delta r=\pm 0.2$. The hard scales vary from $0.75t$ to $1.25t$ to characterize the energy release in decay process. It is necessary to stress that the second uncertainty from $r^{T}(R)$ is absent for the $S$-wave resonance contributions in Table 2. The uncertainties stemming from the weight coefficients $c_{R}$ are not shown explicitly in these Tables, whose effect on the branching ratios via the relation of $\mathcal{B}\propto|c_{R}|^{2}$. For the $S$ and $D$-waves resonance contributions, the twist-3 in the two-kaon DAs are taken as the asymptotic forms for lack of better results from nonperturbative methods, which may give significant uncertainties. We have checked the sensitivity of our results to the choice of the meson radius parameter $r$ [ see Eq. (9)] in the BW parametrization. The variation of its value from 0 $\text{GeV}^{-1}$ to 3.0 $\text{GeV}^{-1}$ results in the change of the branching ratios and polarizations by only a few percents. In general, our results are more sensitive to those hadronic parameters.

Before discussing the results of our calculations in detail, we wish to explain the quoted experimental values that appear in these Tables. The measured branching ratio for each resonant component in the concerned decays are calculated by multiplying its fit fraction and the total three-body decay branching ratio. 333So far, only the $B_{s}\rightarrow J/\psi K^{+}K^{-}$ mode is well measured. Its weighted average branching ratio, given by the Particle Data Group (PDG), is $\mathcal{B}(B_{s}\rightarrow J/\psi K^{+}K^{-})=(7.9\pm 0.7)\times 10^{-4}$ pdg2018 , where the statistical and systematic uncertainties are combined in quadrature. The fit fractions of $P$ and $D$-wave resonances are taken from the most recent LHCb experiment jhep080372017 , which superseded the earlier one from prd87072004 . However, in Ref. jhep080372017 , the $S$-wave component is described in a model-independent pattern, making no assumptions of any $f_{0}$ resonant structures. Therefore, we use the $S$-wave $f_{0}(980)$ and $f_{0}(1370)$ fractions from prd87072004 . Note that the $f_{0}(980)$ fraction is strongly parametrization dependent. For instance, the parameter set by $BABAR$ gives a smaller fit fraction $(4.8\pm 1.0)\%$, while the parameter set by LHCb gives a larger value $(12.0\pm 1.8)\%$ [see Table VI of Ref. prd87072004 ]. Therefore, we quote the central values in a wide range according to the two models rather than a central value plus or minus its statistical and systematic uncertainties for the $f_{0}(980)$ resonance in Table 2. For other charmonium channels, the detailed partial wave analysis for determining various resonance fractions are still missing due to a limited number of events. The quasi-two-body branching ratios can be built from product of two two-body branching ratios when available in the narrow-width limit, namely, $\mathcal{B}(B_{s}\rightarrow XR(\rightarrow K^{+}K^{-}))\approx\mathcal{B}(B_{s}\rightarrow XR)\times\mathcal{B}(R\rightarrow K^{+}K^{-})$. For example, we have used the experimental numbers $\mathcal{B}(B_{s}\rightarrow\eta_{c}\phi(1020))=(5.0\pm 0.9)\times 10^{-4}$ pdg2018 ; jhep070212017 and $\mathcal{B}(\phi(1020)\rightarrow K^{+}K^{-})=(49.2\pm 0.5)\%$ pdg2018 to obtain the experimental branching ratio for $\mathcal{B}(B_{s}\rightarrow\eta_{c}\phi(1020)(\rightarrow K^{+}K^{-}))=(2.5\pm 0.4)\times 10^{-4}$, which is shown in Table 3.

It is clear that the predicted branching ratios of resonant components are consistent with the data except for the tensor resonance $f_{2}(1270)$. From Table. 4, one can see that the predicted branching ratio of $B_{s}\rightarrow J/\psi f_{2}(1270)(\rightarrow K^{+}K^{-})$ is two order of magnitude smaller than the data. We argue that the fit fraction of the $f_{2}(1270)$ component in $B_{s}\rightarrow J/\psi K^{+}K^{-}$ mode jhep080372017 seems to be puzzling since it shows a tension with the corresponding one in $B_{s}\rightarrow J/\psi\pi^{+}\pi^{-}$ 190305530 . 444From discussions with Liming Zhang and Xuesong Liu, the $f_{2}(1270)$ fraction in $B_{s}\rightarrow J/\psi K^{+}K^{-}$ jhep080372017 could be too high because the misidentified background from $B_{d}\rightarrow J/\psi K^{+}\pi^{-}$ in the $f_{2}(1270)$ region may give some systematic uncertainties. For illustration we have explicitly written the relative ratio of $\mathcal{B}(B_{s}\rightarrow J/\psi f_{2}(1270)(\rightarrow K^{+}K^{-}))$ compared to $\mathcal{B}(B_{s}\rightarrow J/\psi f_{2}(1270)(\rightarrow\pi^{+}\pi^{-}))$ in the narrow-width limit as

[TABLE]

in which the common term $\mathcal{B}(B_{s}\rightarrow J/\psi f_{2}(1270))$ in the numerator and denominator cancel out. It is well known that the dominant decay mode of $f_{2}(1270)$ is $\pi\pi$ rather than the $K\bar{K}$, we can thus infer that $\mathcal{R}$ should typically be much less than 1. More specifically, by using the numbers $\mathcal{B}(f_{2}(1270)\rightarrow K\bar{K})=4.6\%$ and $\mathcal{B}(f_{2}(1270)\rightarrow\pi\pi)=84.2\%$ from PDG pdg2018 , we further get $\mathcal{R}=0.04$. Conversely, from the Table 5 in Ref. 190305530 , Table 4 in Ref. prd89092006 , and Table 3 in Ref. jhep080372017 , one can estimate its range from 1.9 to 4.4. It clearly indicates that future improved measurements should take the discrepancy into account. Assuming that the $f_{2}(1270)$ fraction in $B_{s}\rightarrow J/\psi\pi^{+}\pi^{-}$ mode 190305530 is precise enough, combining with above ratio $\mathcal{R}=0.04$, we can infer $\mathcal{B}_{\text{exp}}(B_{s}\rightarrow J/\psi f_{2}(1270)(\rightarrow K^{+}K^{-}))=2.7\times 10^{-7}$. One can see from Table 4 that the predicted branching ratio $\mathcal{B}(B_{s}\rightarrow J/\psi f_{2}(1270)(\rightarrow K^{+}K^{-}))\sim 2.3\times 10^{-7}$ is in agreement with the experiment.

In order to verify the validity of our numerical results, we perform a set of cross-checks.

• (I)

Using our values from Table. 2, we expect that

[TABLE]

where the value of $\mathcal{B}(B_{s}\rightarrow J/\psi f_{0}(980)(\rightarrow\pi^{+}\pi^{-}))=1.15^{+0.52}_{-0.41}\times 10^{-4}$ is read from the previous PQCD calculations prd91094024 . On the experimental side, $BABAR$ measures the ratio of the partial decay width of $f_{0}(980)\rightarrow K^{+}K^{-}$ to $f_{0}(980)\rightarrow\pi^{+}\pi^{-}$ of $0.69\pm 0.32$ using $B\rightarrow KKK$ and $B\rightarrow K\pi\pi$ decays prd74032003 . While BES performs a partial wave analysis of $\chi_{c0}\rightarrow f_{0}(980)f_{0}(980)\rightarrow\pi^{+}\pi^{-}\pi^{+}\pi^{-}$ and $\chi_{c0}\rightarrow f_{0}(980)f_{0}(980)\rightarrow\pi^{+}\pi^{-}K^{+}K^{-}$ in $\psi(2S)\rightarrow\gamma\chi_{c0}$ decay and extracts the ratio as $0.25^{+0.17}_{-0.11}$ prd72092002 . Their weighted average yields $0.35^{+0.15}_{-0.14}$. It can be seen that our estimate in Eq. (36) is consistent with this experimental average value.

• (II)

Combining Tables. 2, 3 and the number in Eq. (36) , we obtain the ratio:

[TABLE]

comply with the latest average of Heavy Flavor Averaging Group (HFAVG) $\mathcal{R}_{f_{0}/\phi}=0.207\pm 0.016$ HFAVG from the measurements plb698115 ; prd85011103 ; plb75684 ; prd84052012

[TABLE]

In comparison to previous theoretical estimation $0.122^{+0.081}_{-0.058}$ obtained in prd83094027 , our value turns out to be larger.

• (III)

Evidence of the $f_{0}(1370)$ resonance in $B_{s}\rightarrow J/\psi\pi^{+}\pi^{-}$ decay is reported by Belle prl106121802 with a significance of $4.2\sigma$. The corresponding product branching fraction is measured to $\mathcal{B}(B_{s}\rightarrow J/\psi f_{0}(1370),f_{0}(1370)\rightarrow\pi^{+}\pi^{-})=3.4^{+1.4}_{-1.5}\times 10^{-5}$ prl106121802 555The PDG also present a value of $\mathcal{B}(B_{s}\rightarrow J/\psi f_{0}(1370)(\rightarrow\pi^{+}\pi^{-}))=4.5^{+0.7}_{-4.0}\times 10^{-5}$ measured by the LHCb Collaboration prd86052006 , which is obtained by multiplying the corresponding normalized fit fraction and the branching ratio of the normalization mode $B_{s}\rightarrow J/\psi\phi(1020)$. Although its central value is consistent with former measurements from Belle, but suffers from sizeable systematic uncertainties. We do not use its result for further calculations.. Combined with our prediction on the $KK$ channel in Table. 2, one can estimate the relative branching ratios of $f_{0}(1370)\rightarrow K^{+}K^{-}/\pi^{+}\pi^{-}$ lie in the range ($0.2\sim 0.5$). Since the situation of the knowledge of the $f_{0}(1370)$ decaying into $KK$ or $\pi\pi$ is rather unclear, above expected values should be investigated further in the future with more precise data.

• (IV)

From Tables 3 and 4, we get another interesting ratio

[TABLE]

in which the two known branching ratios $\mathcal{B}(\phi(1020)\rightarrow K^{+}K^{-})=(49.2\pm 0.5)\%$ and $\mathcal{B}(f^{\prime}_{2}(1525)\rightarrow K^{+}K^{-})=\frac{1}{2}(88.7\pm 2.2)\%$ pdg2018 are used. Our central value is in accordance with the previous theoretical estimation of $0.154^{+0.090}_{-0.070}$ prd95036013 . Experimentally, different Collaborations reported their measurements $\mathcal{R}_{f^{\prime}_{2}/\phi}=0.215\pm 0.049(\text{stat})\pm 0.026(\text{syst})$ (Belle prd88114006 ), $\mathcal{R}_{f^{\prime}_{2}/\phi}=0.264\pm 0.027(\text{stat})\pm 0.024(\text{syst})$ (LHCb prl108151801 ), and $\mathcal{R}_{f^{\prime}_{2}/\phi}=0.19\pm 0.05(\text{stat})\pm 0.04(\text{syst})$ (D0 prd86092011 ). It seems that theoretical predictions are generally smaller than the experimental measurements. None the less, including the errors, both the theoretical predictions and experimental data can still agree with each other.

• (V)

Finally, we estimate the relative branching ratios between two tensor modes

[TABLE]

The current PDG values of

[TABLE]

are dominated by the LHCb measurement prd86052006 ; prd87072004 . Combined with the experiment value $\mathcal{B}(f_{2}(1270)\rightarrow\pi^{+}\pi^{-})=\frac{2}{3}(84.2^{+2.9}_{-0.9})\%$, one obtains the measured ratio $\mathcal{R}_{f_{2}/f^{\prime}_{2}}=0.02\pm 0.01$, which is only half of our prediction in Eq. (39). However, the datum for $f_{2}(1270)$ mode has been further reviewed in Ref 190305530 , the updated branching ratio is $\mathcal{B}(B_{s}\rightarrow J/\psi f_{2}(1270)(\rightarrow\pi^{+}\pi^{-}))=(6.8\pm 1.0)\times 10^{-6}$ with statistical uncertainty only, corresponding to $\mathcal{R}_{f_{2}/f^{\prime}_{2}}=0.05\pm 0.01$. It is clear that our prediction on this ratio is marginally consistent with the updated experiment. In addition, based on the chiral unitary approach for mesons, the authors of Ref. prd90094006 present a larger value $\mathcal{R}_{f_{2}/f^{\prime}_{2}}=0.084\pm 0.046$. Recalling that the theoretical errors are relatively large, so within a $1\sigma$ tolerance, one still can count them as being consistent.

As seen in Table 2, the sum of resonance contributions from $f_{0}(980)$ and $f_{0}(1370)$ is somewhat larger than the $S$-wave total contribution due to the destructive interference between the two resonances. In fact, the best fit model from the LHCb experiment prd87072004 also shows that the destructive interference between $f_{0}(980)$ and $f_{0}(1370)$ resonances in the $B_{s}\rightarrow J/\psi K^{+}K^{-}$ channel. The interference between the two $P$-wave resonances $\phi(1020)$ and $\phi(1680)$ is rather small due to the relatively narrow width of the former ($\Gamma_{\phi(1020)}=$4.25 MeV). Since the contribution of the latter is an order of magnitude smaller, the $P$-wave resonance contribution is almost equal to the $\phi(1020)$ one. By the same token, the $D$-wave resonance contribution mainly come from the $f_{2}^{\prime}(1525)$ component, while other tensor resonance contributions are at least one order smaller. The peak of the high-mass vector resonance $\phi(1680)$ lie almost on the upper limit of the allowed phase space for the $2S$ charmonium modes, their rates suffer a strong suppression and are smaller than that of ground state charmonium channels by almost a factor of 10. Higher-mass $K^{+}K^{-}$ resonances like $f_{2}(1750)$ and $f_{2}(1950)$ are beyond the invariant mass spectra for the $2S$ charmonium modes, their contributions are absent in the last two rows of Table 4. As stated above, any interference contribution between different spin-$J$ states integrates to zero. Therefore, summing over the contributions of the various partial waves, we can obtain the total three-body decay branching ratios

[TABLE]

in which all the uncertainties have been added in quadrature. For the channel $B_{s}\rightarrow J/\psi K^{+}K^{-}$, the obtained branching ratio is slightly smaller than the current PDG average value of $(7.9\pm 0.7)\times 10^{-4}$ pdg2018 . Moreover, keeping in mind that we are not including the nonresonant $S$-wave contribution in our calculations. The small gap might be offset by the nonresonant term and its interference with the resonant components. For other modes, their branching ratios can also reach the order of $10^{-4}$, which is large enough to permit a measurement.

In the literatures, most of the theory studies concentrate on several dominant resonant components. For example, the authors of Ref. prd95036013 considered two dominant $\phi(1020)$ and $f^{\prime}_{2}(1525)$ resonances in the $B_{s}\rightarrow J/\psi K^{+}K^{-}$ decay. The predicted resonance contributions as well as the total three-body decay branching ration are $(5.6\pm 0.7)\times 10^{-4}$, $1.8^{+1.1}_{-0.8}\times 10^{-4}$, 666From discussion with Néstor Quintero, there is a typo for the $f^{\prime}_{2}(1525)$ contribution in the Table IV of prd95036013 , its value should be $1.8$ rather than $0.8$, such that the sum in the last column is $9.3$. and $9.3^{+1.3}_{-1.1}\times 10^{-4}$, respectively. Another earlier paper prd89095026 also discuss the concerned decays in the QCD factorization approach. The three-body branching ratio was obtained by applying Dalitz plot analysis to be $\mathcal{B}(B_{s}\rightarrow J/\psi K^{+}K^{-})=(10.3\pm 0.9)\times 10^{-4}$. In a recent paper 190602489 , the authors have performed phenomenological studies on the $B_{s}\rightarrow J/\psi f_{0}(980)$ decay in the two-body PQCD formalism. With the mixing angle between the $f_{0}(500)$ and $f_{0}(980)$ in the quark-flavor basis adopting as $25^{\circ}$, the calculated branching ratio for the two-body channel $B_{s}\rightarrow J/\psi f_{0}(980)$, was converted into quasi-two-body one as $\mathcal{B}(B_{s}\rightarrow J/\psi f_{0}(980)(\rightarrow K^{+}K^{-}))=4.6^{+2.6}_{-2.3}\times 10^{-5}$. Overall, our results are comparable with these theoretical predictions within the error bars.

The differential branching ratios of the considered decays are plotted on $\omega$ in Fig. 3, in which the green, purple, red, blue, orange, cyan, wine, and black lines show the $f_{0}(980)$, $\phi(1020)$, $f_{2}(1270)$, $f_{0}(1370)$, $f^{\prime}_{2}(1525)$, $\phi(1680)$, $f_{2}(1750)$, and $f_{2}(1950)$ resonance contributions, respectively. To see more clearly all the resonance peaks, especially in the region of the $f_{2}(1270)$ resonance, we draw them in both linear (left panels) and logarithmic (right panels) scales for each decay channel. It is clear that an appreciable peak arising from the $\phi(1020)$ resonance, accompanied by $f_{2}^{\prime}(1525)$. Another three resonance peaks of $f_{0}(980)$, $f_{0}(1370)$, and $\phi(1680)$ have relatively smaller strengths than the $f_{2}^{\prime}(1525)$ one, but their broader widths compensate the integrated strengths over the entire phase space. Therefore, the branching ratios of the four components are predicted to be of a comparable size. Apart from above obvious signal peak, there are two visible structures at about 1750 MeV and 1950 MeV in Fig. 3(a) and 3(c), but not in Fig. 3(e) and 3(g) because the two higher mass regions are beyond the $KK$ invariant mass spectra for the $2S$ charmonium state modes. The contributions of the tensor $f_{2}(1270)$, however, can hardly be seen since its strength is found to be compatible with zero and its peak almost overlap with the tail of those higher mass states like $f_{2}(1750)$ and $f_{2}(1950)$. The obtained distribution for the most of resonance contributions to the $B_{s}\rightarrow J/\psi K^{+}K^{-}$ decay agrees fairly well with the LHCb data shown in Fig. 7 of Ref. jhep080372017 , while other predictions could be tested by future experimental measurements.

Let us now proceed to the polarization fractions which are defined as

[TABLE]

with $\sigma=0,\parallel,\perp$ being the longitudinal, perpendicular, and parallel polarizations, respectively.

The PQCD results for the polarization fractions together with the LHCb data, are listed in Table 5. The sources of the errors in the numerical estimates have the same origin as in the discussion of the branching ratios in Table. 3. For most modes, the transverse polarization fraction $f_{T}=f_{\parallel}+f_{\perp}$ and the longitudinal one are roughly equal. Nevertheless, for the $f_{2}(1950)$ mode, the longitudinal polarization fraction is suppressed to $30\%$ owing to a larger $r^{T}(f_{2}(1950))$ in Eq. (III) enhances its transverse polarization contribution. Even so, the longitudinal polarization fraction is still larger than the experimental value. Of course, taking into account both the theoretical and experimental errors, the deviation is less than $3\sigma$.

For the $P$-wave resonant channels, the parallel polarization fractions are slightly smaller than the corresponding perpendicular one in our calculations, while the LHCb’s data show an opposite behavior for the $\phi(1680)$ mode. As pointed out in Ref. prd98113003 , the relative importance of the parallel and perpendicular polarization amplitudes in the $\rho$ channels are sensitive to the two Gegenbauer moments $a_{2}^{a}$ and $a_{2}^{v}$. The similar situation also exist in this work. Strictly speaking, the Gegenbauer moments in two-hadron DAs are not constants, but depend on the dihadron invariant mass $\omega$. However, the explicit behaviors of those Gegenbauer moments with the $\omega$ are still unknown and the available data are not yet sufficiently precise to control their dependence. Here, we do not consider the $\omega$ dependence and assume the Gegenbauer moments for the resonances with same spin are universal. That is to say it is unlikely to accommodate the measured $B_{s}\rightarrow J/\psi\phi(1020),J/\psi\phi(1680)$ parallel and perpendicular polarization simultaneously with the same set of Gegenbauer moments in PQCD. A further theoretical study of the $\omega$ dependence of the Gegenbauer moments will clarify this issue.

For the $D$-wave mode $B_{s}\rightarrow J/\psi f_{2}(1270)$, compared with the data from the LHCb, our predicted longitudinal polarization is smaller while the two transverse ones are larger [see Table 5]. As stressed before, the $f_{2}(1270)$ fit fraction in the $KK$ mode is unexpected, so its polarizations may have a similar situation. In fact, the best fit model from LHCb prd89092006 ; 190305530 on the $B_{s}\rightarrow J/\psi\pi^{+}\pi^{-}$ decay showing the longitudinal polarization for the $f_{2}(1270)$ component is obviously smaller than the transverse ones. As it is hard to understand why the polarization patterns of $f_{2}(1270)$ resonance decaying into $\pi\pi$ and $KK$ pairs are so different, a refined measurement of the $f_{2}(1270)$ contribution to the $J/\psi K^{+}K^{-}$ mode is urgently needed in order to clarify such issue.

So far, there are several literatures prd89094010 ; prd95036013 ; epjc77610 on the calculation of polarization fractions, focusing more on the $\phi(1020)$ and $f^{\prime}_{2}(1525)$ channels. We found numerically that

[TABLE]

It is clear from Table 5 that our calculations are comparable with theirs within errors. Since the higher mass intermediate states in the concerned decays are still received less attention in both theory and experiment, we wait for future comparison.

IV conclusion

In this paper we carry out an systematic analysis of the $B_{s}$ meson decaying into charmonia and $K^{+}K^{-}$ pair by using the PQCD approach. This type of process is expected to receive dominant contributions from intermediate resonances, such as the vector $\phi(1020)$, tensor $f^{\prime}_{2}(1525)$, and scalar $f_{0}(980)$, thus can be considered as quasi-two-body decays. In addition to the three prominent components mentioned above, some significant excitations in the entire $K^{+}K^{-}$ mass spectrum, which have been well established in the $B_{s}\rightarrow J/\psi K^{+}K^{-}$ decay, are also included. These resonances fall into three partial waves according to their spin, namely, $S$, $P$, and $D$-wave states. Each partial wave contribution is parametrized into the corresponding timelike form factor involved in the two-kaon DAs, which can be described by the coherent sum over resonances sharing the same spin. The $f_{0}(980)$ component is described by a Flatté line shape, while other resonances are modeled by the Breit-Wigner function.

After determining the hadronic parameters involved in the two-kaon DAs by fitting our formalism to the available data, we have calculated each resonance contribution in the processes under consideration. It is found that the largest component is the $\phi(1020)$, followed by $f^{\prime}_{2}(1525)$, with others being almost an order of magnitude smaller. The resultant invariant mass distributions for most resonances in the $B_{s}\rightarrow J/\psi K^{+}K^{-}$ decay show a similar qualitative behavior as the LHCb experiment. Since the interference contributions between any two different spin resonances are zero, summing over various partial wave contributions, we can estimate the total three-body decay branching ratios. The obtained branching ratio of the $B_{s}\rightarrow J/\psi K^{+}K^{-}$ decay is in accordance with available experimental data and numbers from other approaches. The modes involving 2S charmonium have sizable three-body branching ratios, of order $10^{-4}$, which seem to be in the reach of future experiments. As a cross-check, we have discussed some interesting relative branching ratios and compared with available experimental data and other theoretical predictions.

Three polarization contributions were also investigated in detail for the vector-vector and vector-tensor modes. For most of channels, the transverse polarization is found to be of the same size as the longitudinal one and the parallel and perpendicular polarizations are also roughly equal, while for some higher resonance modes, the polarization patterns can be different. The obtained results can be confronted to the experimental data in the future.

Finally, we emphasize that further experimental investigations on the $f_{2}(1270)$ component in the $B_{s}\rightarrow J/\psi K^{+}K^{-}$ decay based on much larger data samples are urgently necessary.

Acknowledgements.

We acknowledge Hsiang-nan Li for helpful discussions and Liming Zhang for enlightening discussions concerning the experiments. This work was supported in part by the National Natural Science Foundation of China under Grants No.11605060 and No.11547020, in part by Natural Science Foundation of Hebei Province under Grant no. A2019209449, and in part by the Program for the Top Young Innovative Talents of Higher Learning Institutions of Hebei Educational Committee under Grant No. BJ2016041. Ya Li is also supported by the Natural Science Foundation of Jiangsu under Grant No. BK20190508.

Appendix A the decay Amplitudes for $B_{s}\rightarrow\eta_{c}R(\rightarrow K^{+}K^{-})$

The decay amplitude can be conventionally written as

[TABLE]

with the Cabibbo-Kobayashi-Maskawa (CKM) matrix elements $V_{ij}$ and the Fermi coupling constant $G_{F}$. $\mathcal{F}(\mathcal{M})$ describes the contributions from the factorizable (nonfactorizable ) diagrams in Fig 1. The superscript $LL$, $LR$, and $SP$ refer to the contributions from $(V-A)\otimes(V-A)$, $(V-A)\otimes(V+A)$, and $(S-P)\otimes(S+P)$ operators, respectively. Performing the standard PQCD calculations, one gets the following expressions:

[TABLE]

[TABLE]

[TABLE]

with color factor $C_{F}=4/3$. $f_{\eta_{c}}$ is the decay constant of the $\eta_{c}$ meson. The expressions for the evolution functions $E$, the hard kernels $h$, and the hard scales $t_{a,b,c,d}$ are referred to the Appendix of Ref. prd91094024 . The forms of $\psi^{v}$ are adopted as our previous works prd90114030 ; epjc75293 . It should be stressed that above factorization formulas are the same for the scalar and vector resonances involved modes except for their different two-kaon DAs. The decay amplitude of the tensor modes should be multiplied by an extra factor $\sqrt{2/3}$, which derives from the different definitions of the polarization vector as aforementioned. In addition, we also consider the vertex corrections to the factorizable diagrams in Fig. 1, whose effects are absorbed into the modified Wilson coefficients as usual bbns1 ; bbns2 ; bbns3 . For the calculation of vertex corrections, one refer to prd65094023 ; prd63074011 for details. The characteristic scale $\Lambda^{5}_{\text{QCD}}=0.225$ GeV at next-to-leading order was used in this work.