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

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.

where these states are shown to be generated from the meson-meson interaction. In Refs [13][14][15], the K + K − S-wave contribution in the φ(1020) resonance region is estimated to be of the order 1 − 10%, in agreement with previous measurements from LHCb [1,16], CDF [17], and ATLAS [18]. The significant S-wave effects may affect measurements of the CP violating phase β s [13,14,19].
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) [20,21]. 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 [22,23]. 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 [24][25][26][27][28][29][30][31][32][33][34]. 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 1 : f 0 (980), f 0 (1370), φ(1020), φ(1680), f 2 (1270), f ′ 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 [1,7]. All resonances are commonly described by Breit Wigner (BW) distributions, except for the f 0 (980) state, which is modelled by a Flatté function [35].
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. 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 = M where ζ is the kaon momentum fraction. The mass ratio r = m/M with the charmonium mass m. The factor η = p 2 /(M 2 − m 2 ) is defined in terms of the invariant mass squared of the kaon pair p 2 = ω 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

II. KINEMATICS AND THE TWO-KAON DISTRIBUTION AMPLITUDES
in which x B , x 3 , z denote the longitudinal momentum fractions, and k iT represent the transverse momenta.
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 [27][28][29][30]36], 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. [36][37][38]. Below, we briefly describe the two-kaon DAs in three partial waves and the associated form factors.

A. S-wave two-kaon DAs
The S-wave two-kaon DAs are introduced in analogy with the case of two-pion ones [22,24], which are organized into with the null vectors n = (1, 0, 0 T ) and v = (0, 1, 0 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 [39], we adopt their asymptotic models as shown below [24,30]: with the isoscalar scalar form factor F S (ω 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 [40]. 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 [1]. For the former, we use a Flatté description, while the latter is modelled by BW functions. The scalar form factor F S (ω 2 ) can be written as 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 ρ ππ = 2q π /ω, ρ KK = 2q K /ω, where q π(K) is the pion (kaon) momentum in the dipion (dikaon) rest frame. The exponential factor F KK = e −αq 2 K with α = 2.0 ± 0.25 GeV −2 [41,42] is introduced above the KK threshold and serves to reduce the ρ KK factor as the invariant mass increases. The constants g ππ and g KK are the f 0 (980) couplings to ππ and KK final states respectively. We use g ππ = 167 MeV and g KK /g ππ = 3.47 as determined by LHCb [43]. The BW amplitude in generic form is where m R is the resonance pole mass and Γ(ω 2 ) is its energy-dependent width which may be parametrized in a form that ensures the correct behavior near threshold, Here Γ 0 and q K0 are Γ(ω 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 [44] for scalar, vector and tensor states are with z = r 2 q 2 K and z 0 represents the value of z when ω = m R . The meson radius parameters r are set to 1.5 GeV −1 consistent with the literature [1].

B. P -wave two-kaon DAs
In Ref. [29] 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 where the superscripts L and T on the left-hand side denote the longitudinal polarization and transverse polarization, respectively. Here ǫ µνρσ is the totally antisymmetric unit Levi-Civita tensor with the convention ǫ 0123 = 1. The transverse polarization vector ǫ T for the dikaon system has the same form as that of dipion [29]. The various twists DAs in Eq. (9) can be expanded in terms of the Gegenbauer polynomials: where the two P -wave form factors F P and F ⊥ P 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 [45]. The Gegenbauer moments a i 2 will be regarded as free parameters and determined in this work. As mentioned in the Introduction, the form factor F P is given by the coherence summation of the two vector resonances φ(1020) and φ(1680), For the F ⊥ P , we employ the approximate relation F ⊥ being the tensor (vector) decay constant for the corresponding vector meson.

C. 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 ǫ • which is related to the polarization tensor ǫ µν (λ) with helicity λ [46]. Following a similar procedure, we decompose the D-wave two-kaon DAs associated with longitudinal and transverse polarizations into respectively, where the prefactor 2 3 ( 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 Φ 0 D and Φ T D have similar asymptotic forms as the corresponding ones for a tensor meson. More precisely, Here, we solely employ the first nonvanishing leading term in the expansion for previously mentioned reasons. The moments a 0 1 and a T 1 are regarded as free parameters and determined in the next section. Note that the ζ dependent terms L(ζ) and T (ζ) are different from those in Eq. (10). We will derive their expressions later. The kaon tensor form factor F D (ω 2 ) can be represented by where the summation is performed over the intermediate tensor mesons: f ′ 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 [47,48] Next we derive the ζ 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 [46,47] where p and m T are momentum and mass of the tensor meson, respectively. The two interpolating currents j µν (0) and j µνρ (0) are defined in [46,47]. Let us begin with the local matrix element K + (p 1 )K − (p 2 )|j µν (0)/j µνρ (0)|0 associated with the D-wave form factors. Under the tensor-meson-dominant hypothesis [47], inserting the tensor intermediate in above matrix element, we get with D f2 the resonance propagator [49]. The coupling constant g f2KK is defined by the matrix element [47]. When applying the formula Eq. (16) and the completeness relation Eq. (17) then leads to the following equations explicitly where power suppressed terms of order m 2 T in the parenthesis have been omitted. Utilizing the approximation relation in which the coefficient 1 − 6ζ + 6ζ 2 is absorbed into the longitudinal polarization DAs, giving rise to its ζ dependence. The matrix element in Eq. (19) for the choice µ, ν, ρ = +, −, x is proportional to 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 ζ dependent factors of the longitudinal and transverse polarization DAs can be written as respectively. Above expressions can also be checked from the partial wave expansion of the helicity amplitudes. As is well known, the helicity 0 component is expanded in terms of Legendre polynomials P l (cosθ), while the helicity ±1 ones proceeding in derivatives of the Legendre polynomials P ′ l (cosθ) [15]. For l = 2 D-wave amplitudes, the helicity angle θ is encoded into the Wigner-d functions, schematically [2]: Following a similar prescription in the dipion system [50,51], the polar angle θ of the K + in the rest frame of the dikaon is related to ζ through the relation with m K the kaon mass. Neglecting the kaon mass and employing Eq. (23), we also arrive at Eq. (22).

D. The differential branching ratio
The double differential branching ratio reads [2]  Belle [54] where the differential variable dcosθ is replaced by dζ via Eq. (24) in the limit of massless. The three-momenta of the kaon and charmonium in the rest reference frame of the KK system are given by respectively, with the standard Källén function λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + ac + bc). The complete amplitude A through resonance intermediate for the concerned decay is decomposed into where A S , A P , and A D denote the corresponding three partial wave decay amplitudes. Since the ζ-dependent terms appear as an overall factor in each partial wave decay amplitudes, integrating the double differential distribution of Eq. (25) over ζ gives for the differential invariant mass branching ratio where the factors 1/3, 1/5, · · · are extracted from the individual helicity amplitudes for the integral of ζ. The terms A 0 , A , and A ⊥ represent the longitudinal, parallel, and perpendicular polarization amplitudes in the transversity basis, respectively 2 . Note that interference between different partial wave vanishes because the ζ functions in Eqs. (5), (10), and (22), corresponding to S, P , and D partial waves, are orthogonal. In the J/ψ and ψ(2S) cases, the decay amplitudes A S and A P here can be straightforwardly obtained from the previous publications [24,29] by replacing the two-pion form factors and all pion masses and momenta with the respective kaon quantities. For A D , its factorization formula can be related to A P by making the following replacement, For the involved η c and η 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) [2]: II: Branching ratios of S-wave resonance contributions to the Bs → (J/ψ, ψ(2S), ηc, ηc(2S))K + K − decays. Theoretical errors correspond to the uncertainties of Gegenbauer moments and hard scales, respectively.
a We are not including contributions from the nonresonant S-wave. b We quote the range of measurement since the fit fraction of f 0 (980) is strongly parametrization dependent. c The fit fraction statistical and systematic are added in quadrature.
The information on the decay constants (in units of GeV), the Wolfenstein parameters, together with the lifetime of B s mesons are adopted as [29,[36][37][38]47] f The masses and widths of the BW resonances are listed in Table I, while the Flatté parameters for the f 0 (980) have been given in the previous section.
As mentioned before, the form factors ratio r T (R) = F ⊥ /F is approximately equal to the ratio of two decay constants f T /f . From the numbers in Eq. (31), we have 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 factors c i , Gegenbauer moments a i , and form factor ratios r T (R), in each partial wave. Their central values are fixed to be S-wave : c f0(1370) = 0.12e i π 2 , a 1 = 0.8, P-wave : where the phases of c f2(1270) and c f2(1750) are set to π to ensure the normalization of the form factor F D (0) = 1. When fitting to experimental data, we assume that the Gegenbauer moments are the same for all the resonances with the same spin. We first take into account three largest components in their individual partial wave contribution to fix the corresponding Gegenbauer moments, and then determine r T (R) and the module of c R for each resonance according to its polarization fractions and fit fraction from the LHCb [7], respectively. The calculated branching ratios of S, P , and D-wave resonance contributions to the B s → (J/ψ, ψ(2S), η c , η c (2S))K + K − decays are collected in Tables II, III, (32) and (33) are general assigned to be δr = ±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 II. 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 B ∝ |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. 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 III: Branching ratios of P -wave resonance contributions to the Bs → (J/ψ, ψ(2S), ηc, η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 [2,7,55] are combined in quadrature.

Modes
B(R = φ(1020)) B(R = φ(1680)) P-wave The fit fractions determined from the Dalitz plot analysis have been converted into the branching ratio measurements. b The experimental data is obtained by the product of B(Bs → ηcφ(1020)) and B(φ(1020) → K + K − ). c The experimental data is obtained by the product of B(Bs → ψ(2S)φ(1020)) and B(φ(1020) → K + K − ).  Table. III. The statistical and systematic uncertainties from [7] are combined in quadrature.
calculated by multiplying its fit fraction and the total three-body decay branching ratio 3 . The fit fractions of P and D-wave resonances are taken from the most recent LHCb experiment [7], which superseded the earlier one from [1]. However, in Ref. [7], 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 [1]. 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 ± 1.0)%, while the parameter set by LHCb gives a larger value (12.0 ± 1.8)% [see Table VI of Ref. [1]]. 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 II. 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, B(B s → XR(→ K + K − )) ≈ B(B s → XR) × B(R → K + K − ). For example, we have used the experimental numbers B(B s → η c φ(1020)) = (5.0 ± 0.9) × 10 −4 [2,55] and B(φ(1020) → K + K − ) = (49.2 ± 0.5)% [2] to obtain the experimental branching ratio for B(B s → η c φ(1020)(→ K + K − )) = (2.5 ± 0.4) × 10 −4 , which is shown in Table III. 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. IV, one can see that the predicted branching ratio of B s → J/ψf 2 (1270)(→ K + K − ) is two order of magnitude smaller than the data. We argue that the fit fraction of the f 2 (1270) component from Ref. [7] seems to be questionable 4 . It is worth noting that there are some tension in the LHCb fitting between B s → J/ψπ + π − [56] and B s → J/ψK + K − [7]. For illustration we have explicitly written the relative ratio of B(B s → J/ψf 2 (1270)(→ K + K − )) compared to B(B s → J/ψf 2 (1270)(→ π + π − )) in the narrow-width limit as in which the common term B(B s → J/ψf 2 (1270)) in the numerator and denominator cancel out. It is well known that the dominant decay mode of f 2 (1270) is ππ rather than the KK, we can thus infer that R should typically be much less than 1. More specifically, by using the numbers B(f 2 (1270) → KK) = 4.6% and B(f 2 (1270) → ππ) = 84.2% from PDG [2], we further get R = 0.04. Conversely, from the Table 5 in Ref. [56], Table 4 in Ref. [42], and Table 3 in Ref. [7], 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 → J/ψπ + π − mode [56] is precise enough, combining with above ratio R = 0.04, we can infer B exp (B s → J/ψf 2 (1270)(→ K + K − )) = 2.7 × 10 −7 , which is consistent with our prediction in Table IV.
In order to verify the validity of our numerical results, we perform a set of cross-checks.
In comparison to previous theoretical estimation 0.122 +0.081 −0.058 obtained in [64], our value turns out to be larger. (III) Evidence of the f 0 (1370) resonance in B s → J/ψπ + π − decay is reported by Belle [65] with a significance of 4.2σ. The corresponding product branching fraction is measured to B(B s → J/ψf 0 (1370), f 0 (1370) → π + π − ) = 3.4 +1.4 −1.5 ×10 −5 [65] 5 . Combined with our prediction on the KK channel in Table. II, one can estimate the relative branching ratios of f 0 (1370) → K + K − /π + π − lie in the range (0.2 ∼ 0.5). Since the situation of the knowledge of the f 0 (1370) decaying into KK or ππ is rather unclear, above expected values should be investigated further in the future with more precise data.
As seen in Table II, 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 [1] also shows that the destructive interference between f 0 (980) and f 0 (1370) resonances in the B s → J/ψK + K − channel. The interference between the two P -wave resonances φ(1020) and φ(1680) is rather small due to the relatively narrow width of the former (Γ φ(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 φ(1020) one. By the same token, the D-wave resonance contribution mainly come from the f ′ 2 (1525) component, while other tensor resonance contributions are at least one order smaller. The peak of the high-mass vector resonance φ(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 IV. 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 in which all the uncertainties have been added in quadrature. For the channel B s → J/ψK + K − , the current PDG average value is (7.9 ± 0.7) × 10 −4 [2], which is slightly larger than our calculation in Eq. (40). 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. [9] considered two dominant φ(1020) and f ′ 2 (1525) resonances in the B s → J/ψK + K − decay. The predicted resonance contributions as well as the total three-body decay branching ration are (5.6 ± 0.7) × 10 −4 , 1.8 +1.1 −0.8 × 10 −4 6 , and 9.3 +1.3 −1.1 × 10 −4 , respectively, which is comparable to our calculations. Another earlier paper [8] also discuss the concerned decays in the QCD factorization approach. The three-body branching ratio was obtained by applying Dalitz plot analysis to be B(B s → J/ψK + K − ) = (10.3 ± 0.9) × 10 −4 , which is a little larger than our prediction. In a recent paper [66], the authors have performed phenomenological studies on the B s → J/ψf 0 (980) decay in the two-body PQCD formalism. The calculated branching ratio for the two-body channel B s → J/ψf 0 (980) was converted into quasi-two-body one as B(B s → J/ψf 0 (980)(→ K + K − )) = 4.6 +2.6 −2.3 × 10 −5 , which agree well with our prediction. The overall agreement between two-body and three-body frameworks implies that the PQCD approach is a consistent theory for the hadronic charmonium B s decays.  V: Polarization fractions for the decays Bs → (J/ψ, ψ(2S))φ/f2(→ K + K − ). For theoretical errors, see Table. III. The experimental data are taken from Ref. [7], where the statistical and systematic uncertainties are combined in quadrature.
Modes The differential branching ratios of the considered decays are plotted on ω in Fig. 2, in which the green, purple, red, blue, orange, cyan, wine, and black lines show the f 0 (980), φ(1020), f 2 (1270), f 0 (1370), f ′ 2 (1525), φ(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 φ(1020) resonance, accompanied by f ′ 2 (1525). Another three resonance peaks of f 0 (980), f 0 (1370), and φ(1680) have relatively smaller strength than the f ′ 2 (1525) one, but their broader widths compensate the integrated strength 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. 2(a) and 2(c), but not in Fig.  2(e) and 2(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 → J/ψK + K − decay agrees fairly well with the LHCb data shown in Fig. 7 of Ref. [7], while other predictions could be tested by future experimental measurements.
Let us now proceed to the polarization fractions which are defined as with σ = 0, , ⊥ being the longitudinal, perpendicular, and parallel polarizations, respectively. The PQCD results for the polarization fractions together with the LHCb data, are listed in Table V. The sources of the errors in the numerical estimates have the same origin as in the discussion of the branching ratios in Table. III. For most modes, the transverse polarization fraction f T = f +f ⊥ 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. (33) 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σ.
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 φ(1680) mode. As pointed out in Ref. [29], the relative importance of the parallel and perpendicular polarization amplitudes in the ρ channels are sensitive to the two Gegenbauer moments a a 2 and a v 2 . 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 ω. However, the explicit behaviors of those Gegenbauer moments with the ω are still unknown and the available data are not yet sufficiently precise to control their dependence. Here, we do not consider the ω 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 → J/ψφ(1020), J/ψφ(1680) parallel and perpendicular polarization simultaneously with the same set of Gegenbauer moments in PQCD. A further theoretical study of the ω dependence of the Gegenbauer moments will clarify this issue.
For the D-wave mode B s → J/ψ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 V]. 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 [42,56] on the B s → J/ψπ + π − 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 ππ and KK pairs are so different, a refined measurement of the f 2 (1270) contribution to the J/ψK + K − mode is urgently needed in order to clarify such issue.
So far, there are several literatures [9,36,68] on the calculation of polarization fractions, focusing more on the φ(1020) and f ′ 2 (1525) channels. We found numerically that J/ψφ(1020) : which are in agreement with ours as well as the experimental data. 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 important contributions from intermediate resonances, such as the vector φ(1020), tensor f ′ 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 → J/ψ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 with 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 φ(1020), followed by f ′ 2 (1525), with others being almost an order of magnitude smaller. The resultant invariant mass distributions for most resonances in the B s → J/ψ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 get the total three-body decay branching ratios. The obtained branching ratio of the B s → J/ψ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.