Resonant contributions to three-body $B_{(s)} \to [ D^{(*)}, \bar{D}^{(*)} ] K^+K^-$ decays in the perturbative QCD approach

In this work, we study the $S$, $P$ and $D$ wave resonance contributions to three-body decays $B_{(s)} \to [ D^{(*)}, \bar{D}^{(*)} ] K^+K^-$ by employing the perturbative QCD (PQCD) approach, where the kaon-kaon invariant mass spectra are dominated by the $f_0(980),f_0(1370),$ $\phi(1020),\phi(1680),f_2(1270),f^{\prime}_2(1525),f_2(1750)$ and $f_2(1950)$ resonances. The $KK$ $S$-wave component $f_0(980)$ is modeled with the Flatt\'e formalism, while other resonances are described by the relativistic Breit-Wigner (BW) line shape. The corresponding decay channels are studied by constructing the kaon-kaon distribution amplitude $\Phi_{KK}$, which captures important final state interactions between the kaon pair in the resonant region. We found that the PQCD predictions for the branching ratios for most considered decays agree with currently available data within errors. The associated polarization fractions of those vector-vector and vector-tensor decay modes are also predicted, which are expected to be tested in the near future experiments. The invariant mass spectra for the corresponding resonances in the $B_{(s)} \to [D^{(*)}, \bar{D}^{(*)} ] K^+K^-$ decays are well established, which can be confronted with the precise data from the LHCb and Belle II experiments.

On the theoretical side, the two-body charmed decays B (s) → [D ( * ) ,D ( * ) ][S, P, V, T ] (here S, P, V and T denote the scalar, pseudoscalar, vector and tensor mesons) have been investigated within the framework of the PQCD factorization approach [20][21][22][23][24][25][26]. The channels induced by the b → c transitions are CKM favored, while those induced by the b → u transitions are CKM suppressed. Thus, the b → u decays will have smaller branching ratios. The interference between the b → c and b → u transitions gives the measurement of the CKM angle γ. As is well-known, the charmed decays of B (s) are more complicated because of the hierarchy of the scale involved compared with the decays of B (s) mesons to the light vector mesons. For example, the B → D transitions involve three scales: the B meson mass m B , the D meson mass m D , and the mass difference from the heavy meson and the heavy quarkΛ = m B − m b ∼ m D − m c , which are strikingly different from each other. Although, the factorization has been proved in soft-collinear effective theory [27], it needs more inputs than the PQCD approach. It can be found that the momentum square of the hard gluon connecting the spectator quark is only a factor of (1 − m 2 D /m 2 B ) to that of the B → light transitions for B → D transitions, which ensures that PQCD can also work well in B → D transitions. There are also many other traditional methods and approaches to estimate the B → D transitions, such as the heavy quark effective theory (HQET) [28,29], light cone sum rules (LCSR) [30][31][32] and lattice QCD (LQCD) [33][34][35].
As addressed before, the B → DKK decay is expected to proceed through φ → KK intermediate state. Moreover, this process can also be dominated by a series of other resonances in S, P , and D waves. In this work, we will study the S, P and D wave resonance contributions to three-body decays B → DKK by employing the PQCD approach, where the kaon-kaon invariant mass spectra are dominated by the f 0 (980), f 0 (1370), φ(1020), φ(1680), f 2 (1270), f ′ 2 (1525), f 2 (1750) and f 2 (1950) resonances. However, the theory of three-body non-leptonic decays is still in an early stage of development. Three-body B meson decay modes do receive the entangled resonant and nonresonant contributions, as well as the possible finalstate interactions (FSIs) [36][37][38], whereas the relative strength of these contributions vary significantly in different regions of the Dalitz plots [39,40]. In this respect, three-body decays are considerably more challenging than two-body decays, but provide a number of theoretical and phenomenological advantages. On the one hand, the number of different three-body final states is about ten times larger than the number of two-body decays. On the other hand, each final state has a non-trivial kinematic multiplicity (a twodimensional phase space) as opposed to two-body decays where the kinematics is fixed by the masses. This leads to a much richer phenomenology, but there is no proof of factorization for the three-body B decays at present. We can only restrict ourselves to specific kinematical configurations on basis of the Dalitz plot analysis. The Dalitz plot contains different regions with "specifical" kinematics [41,42]. The central region corresponds to the case where all three final particles fly apart with similar large energy and none of them moves collinearly to any other. This situation contains two hard gluons and is power and α s suppressed with respect to the amplitude at the edge. The corners correspond to the case in which one final particle is approximately at rest (i.e. soft), and the other two are fast and back-to-back. The central part of the edges corresponds to the case in which two particles move collinearly and the other particle recoils back: this is called the quasi-two-body decay. This situation exists particularly in the low ππ or Kπ or KK invariant mass region of the Dalitz plot. Thereby, it is reasonable to assume the validity of factorization for the quasitwo-body B meson decay. Naturally the dynamics associated with the pair of final state mesons can be factorized into a two-meson distribution amplitude (DA) Φ h 1 h 2 [43][44][45][46][47][48][49].
In recent years, based on the PQCD approach, more and more detailed analysis on the three-body B meson hadronic decays have been performed in the low energy resonances on ππ, KK , Kπ and πη channels [50][51][52][53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68]. Other theoretical approaches for describing the three-body hadronic decays of B mesons based on the symmetry principles and factorization theorems have been developed. The QCD-improved factorization (QCDF) [69][70][71][72] has been widely adopted in the study of the three-body charmless hadronic B meson decays [42,[73][74][75][76][77][78][79][80][81][82][83][84]. The U -spin and flavor SU (3) symmetries were used in Refs. [85][86][87][88][89][90]. Unlike the collinear factorization in the QCD factorization approach and soft-collinear effective theory, the k T factorization is utilized in the PQCD approach. In this approach, the transverse momentum of valence quarks in the mesons is kept to avoid the endpoint singularity [91,92]. The Sudakov factors from the k T resummation have been included to suppress the long-distance contributions from the large-b region with b being a variable conjugate to k T . Therefore, one can calculate the color-suppressed channels as well as the color-allowed channels in charmed B decays within the PQCD approach. The conventional non-calculable annihilation-type decays are also calculable in the PQCD approach, which is proved to be the dominant strong phase in B decays for the direct CP asymmetry.
As aforementioned for the cases of the quasi-two-body decays, the two mesons (h 1 h 2 ) move collinearly fast, and the bachelor meson h 3 is also energetic and recoils against the meson pair in the B meson rest frame in the quasi-two-body B → (h 1 h 2 )h 3 decays. The interaction between the meson pair and the bachelor meson is regarded as to be power suppressed. The typical PQCD factorization formula for the B → (h 1 h 2 )h 3 decay amplitude can be described as the form of [50], where H is the hard kernel, Φ B and Φ h 3 are the universal wave functions of the B meson and the bachelor meson, respectively. The hard kernel H describes the dynamics of the strong and electroweak interactions in three-body hadronic decays in a similar way as the one for the corresponding two-body decays. The wave functions Φ B and Φ h 3 absorb the non-perturbative dynamics in the process. The Φ h 1 h 2 is the two-hadron (KK pair in this work) DA, which involves the resonant and nonresonant interactions between the two moving collinearly mesons. The present paper is organized as follows. In Sec. II, we give a brief introduction for the theoretical framework and the total decay amplitudes with the wilson coefficients, CKM matrix elements and the amplitudes of four-quark operators needed in the calculation will be given in Sec. III. Section IV contains the numerical values and some discussions. A brief summary is given in section V. The Appendix collects the explicit PQCD factorization formulas for all the decay amplitudes.

A. The effective Hamiltonian and kinematics
For B → D ( * ) (R →)KK decays, the weak effective Hamiltonian can be specified as follows [93]: where V cb(q) and V ub(q) are the CKM matrix elements and R denotes the various partial wave resonances f 0 (980), f 0 (1370), φ(1020), φ(1680), f 2 (1270), f ′ 2 (1525), f 2 (1750) and f 2 (1950) 1 , respectively. The explicit expressions of the local four-quark tree operators O 1,2 (µ) and the corresponding Wilson coefficients C 1,2 (µ) can be found in Ref. [93]. The q in Eq. (2) represents the quark d or s. Noted that only tree diagrams contribute to these processes, which shows that there is no direct CP asymmetry in these decays. The typical Feynman diagrams for the decays B (s) →D ( * ) (R →)KK and B (s) → D ( * ) (R →)KK are shown in Fig. 1 and 2, respectively.
We will work in the B meson rest frame and employ the light-cone coordinates for momentum variables. In the light-cone coordinates, we let the kaon pair and the final-state D ( * ) move along the directions n = (1, 0, 0 T ) and v = (0, 1, 0 T ), respectively. The B meson momentum p B , the total momentum of the kaon pair, p = p 1 + p 2 , the final-state D ( * ) momentum p 3 , and the quark momentum k i in each meson are defined in the following form where m B is the mass of the B meson, η = is the mass of the bachelor meson, and the invariant mass squared ω 2 = (p 1 + p 2 ) 2 = p 2 . The momentum fractions x B , z, and x 3 run from zero to unity, respectively. In the heavy quark limit, the mass differenceΛ between b-quark (c-quark) and B(D) meson is negligible.
As usual we also define the momentum p 1 and p 2 of kaon pair as with ζ = p + 1 /P + characterizing the distribution of the longitudinal momentum of the kaon and B. Wave functions of B meson and the D ( * ) mesons The light-cone matrix element of the B meson can be decomposed as [94]: where q represents u or d or s quark. According to the above equation, there are two different wave functions φ B andφ B in the B meson distribution amplitudes, which obey the following normalization conditions: In general, one should consider these two Lorentz structures in the calculations of B decays as shown in Refs. [95]. However, we neglect the contribution ofφ B because of the numerical suppression in this work. Then, the wave function of B meson can be written as [94,[96][97][98][99][100][101] with the widely used B-meson DA in the PQCD approach [96,98] where the normalization factor N B depends on the values of ω B and f B and defined through the normalization relation . ω B is a free parameter and ω B = 0.40 ± 0.04 GeV and ω Bs = 0.50 ± 0.05 GeV [96,100,101] are used in the numerical calculations. Very recently, a new method was proposed to calculate the B meson light-cone DA from lattice QCD, which can be used as an updated input for the B meson DA in the future [102].
For the D ( * ) meson, in the heavy quark limit, the two-parton light-cone DA can be written as [20][21][22][23][24][25] where with C D = 0.5 ± 0.1, ω = 0.1 GeV and f D + = 211.9 MeV [103] for D meson. In the above models, x is the momentum fraction of the light quark in D meson. We determine the decay constant of the vector meson D * by using the relation f D * = m D m D * f D based on the heavy quark effective theory [104].

C. Two-kaon DAs
Below, we briefly introduce the S, P , and D-wave two-kaon DAs and the corresponding time-like form factors used in our framework. It will be shown that resonant contributions through two-body channels can be included by parameterizing the two-kaon DAs. The S-wave two-kaon DAs are described in the following form [61], In what follows the subscripts S, P , and D are always associated 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 [105]. For the scalar resonances f 0 (980) and f 0 (1370), the asymptotic forms of the individual DAs in Eq. (12) have been parameterized as [43][44][45][46] with the time-like scalar form factor F S (ω 2 ) and the Gegenbauer coefficient a S . The P -wave two-pion DAs related to both longitudinal and transverse polarizations have been studied in Ref. [106]. Naively, the P -wave two-kaon ones can be obtained by replacing the pion vector form factors with the corresponding kaon ones. The explicit expressions of the P -wave kaon-kaon DAs associated with longitudinal (L) and transverse (T) polarization are described as follows, The two-kaon DAs for various twists are expanded in terms of the Gegenbauer polynomials: with the two P -wave form factors F P (ω 2 ) and F ⊥ P (ω 2 ) and the Gegenbauer coefficients a i 2P . We introduce the D-wave two-kaon DAs associated with longitudinal and transverse polarizations as follows [61], The D-wave DAs are given as with the ζ dependent factor P 2 (2ζ − 1) = 1 − 6ζ(1 − ζ) and T (ζ) = (2ζ − 1) ζ(1 − ζ). F ,⊥ D (ω 2 ) are the D-wave time-like form factors and the Gegenbauer moments a 0,T D have been determined in our previous work [61].
The strong interactions between the resonance and the final-state meson pair, including elastic rescattering of the final-state meson pair, can be factorized into the time-like form factor F S,P,D (ω 2 ), which is guaranteed by the Watson theorem [107]. For a narrow resonance, we usually use the relativistic BW line shape to parameterize the time-like form factor F (ω 2 ). The explicit expression is [108], where the corresponding weight coefficients c i are determined based on the normalization condition F (0) = 1. The m i and Γ i are the pole mass and width of the corresponding resonances shown in Table I, respectively. The mass-dependent width Γ i (ω) is defined as The | p 1 | is the momentum vector of the resonance decay product measured in the resonance rest frame, while | p 0 | is the value of | p 1 | at ω = m i . The explicit expression of kinematic variables | p 1 | is  [111] with the Källén function λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + ac + bc). L R is the orbital angular momentum in the dikaon system and L R = 0, 1, 2, ... corresponds to the S, P, D, ... partial-wave resonances. Due to the limited studies on the form factor F ⊥ (ω 2 ), we use the two decay constants f (T ) i of the intermediate particle to estimate the form factor ratio r The BW formula does not work well for f 0 (980), because its pole mass is close to the KK threshold. The resulting line shape above and below the threshold of the intermediate particle is called Flatté parametrization [112]. If the coupling of a resonance to the channel opening nearby is very strong, the Flatté parametrization shows a scaling invariance and does not allow for an extraction of individual partial decay widths. If the scalar resonance lies under the KK threshold, the position of the peak in the mass spectrum does not coincide with the pole mass of the resonance. To solve the problem, the finite width to the propagator of scalar resonance has been taken into consideration by Achasov's [113][114][115], which is the one-loop contribution to the self-energy of the scalar resonance from the two-particle intermediate states. More details can be found in Refs. [113][114][115][116][117][118]. In addition, the exponential factor F KK = e −αq 2 K with α = (2.0 ± 0.25)GeV −2 is introduced above the KK threshold and serves to reduce the ρ KK factor as the invariant mass increases, where q K is the koan momentum in the kaon-kaon rest frame [119]. This parametrization decreases the f 0 (980) width above KK threshold slightly. In this work, the invariant mass of the dikaon is above the K + K − threshold, we have tested the impact of the self-energy correction and found it is small. Thus, we employ the modified Flatté model suggested by Bugg [119] following the LHCb collaboration [120,121], The coupling constants g ππ = 0.167 GeV and g KK = 3.47g ππ [120,121] describe the f 0 (980) decay into the final states π + π − and K + K − , respectively. The phase space factors ρ ππ and ρ KK read as [112,120,122] D. The differential branching ratio The double differential branching ratio can be obtained as [103] The three-momenta of the kaon and D meson in the KK center-of-mass frame are given by .
The complete amplitude A through intermediate resonances for the concerned decay channels can be written as the summation of A S , A P and A D : where A S , A P , and A D denote the corresponding three S, P and D wave decay amplitudes. Due to the angular momentum conservation requirement, the vector mesons φ,D * , D * and tensor meson f 2 →]KK decays, both the longitudinal polarization and the transverse polarization contribute. The amplitudes can be decomposed as follows: where A L is the longitudinally polarized decay amplitude, A N and A T are the transversely polarized contributions. Therefore, the total decay amplitude for 2 →]KK decays can be expressed as where A 0 , A and A ⊥ are defined as: The polarization fractions f λ with λ = 0, , and ⊥ are described as with the normalisation relation f 0 + f + f ⊥ = 1.

III. CALCULATION OF DECAY AMPLITUDES IN PQCD APPROACH
In this section, we intend to calculate the relevant decay amplitudes. For the considered B (s) → D ( * ) (R →)KK decays, the analytic formulas for the corresponding decay amplitudes are of the following form: • S-wave: • P-wave: • D-wave: While for B (s) → D ( * ) (R →)KK channels, the total decay amplitudes are written as: • S-wave: • P-wave: • D-wave: in which the combinations a i of the Wilson coefficients are defined as: 2 →]KK decays, the superscript i = L, N, T , represent longitudinal, parallel, and transverse polarization contributions respectively. F LL e(a) and M LL e(a) describe the contributions from the factorizable emission (annihilation) and non-factorizable emission (annihilation) diagrams as show in Fig. 1 and 2. The explicit expressions of the amplitudes F LL e(a) and M LL e(a) will be given in the Appendix.

IV. NUMERICAL RESULTS
In this section, let us firstly list the parameters used in our numerical calculations, such as the masses (in units of GeV) [103]: The values of the Wolfenstein parameters are adopted as given in the Ref. [103]: The decay constants (in units of GeV) and the B meson lifetimes (in units of ps) are chosen as [56,123,124] The form factor ratio r T , the coefficients c i and the Gegenbauer moments a S,P,D are adopted the same values as those determined in Ref. [61] r T (φ(1020)) = 0.865, r T (φ(1680)) = 0.6, r T (f 2 (1270)) = 1.15, By using the Eqs. (37)(38)(39)(40)(41)(42)(43), the decay amplitudes in Sec. III and Appendix and all the input quantities, the resultant branching ratios B and the polarization fractions f λ together with the available experimental measurements for the considered quasi-two-body decays B (s) → [D ( * ) ,D ( * ) ](R →)KK are summarized in Tables II-VI. In our numerical calculations for the branching ratios and polarization fractions, the first theoretical uncertainty results from the parameters of the wave functions of the initial and final states, such as the shape parameter ω B = 0.40 ± 0.04 GeV or ω Bs = 0.50 ± 0.05 GeV and C D = 0.5 ± 0.1, and the decay constant f B (s) in B (s) meson wave function. The second one is due to the Gegenbauer moments in DAs of KK pair with different intermediate resonances, which are supposed to be varied with a 20% range for the error estimation. With the improvements of the experiments and the deeper theoretical developments, this kind of uncertainties will be reduced. The last one is caused by the variation of the hard scale t from 0.75t to 1.25t (without changing 1/b i ) and the QCD scale Λ QCD = 0.25 ± 0.05 GeV, which characterizes the effect of the next-to-leading-order QCD contributions. The possible errors due to the uncertainties of m c and CKM matrix elements are very small and can be neglected safely.
Compared with B (s) → D ( * ) R → D ( * ) KK decays, the B (s) →D ( * ) R →D ( * ) KK ones are enhanced by the CKM matrix elements |V cb /V ub | 2 , especially for those without a strange quark in the four-quark operators. So for most of the B (s) →D ( * ) R →D ( * ) KK decays, the branching ratios are at the order from 10 −7 to 10 −5 ; while for the B (s) → D ( * ) R → D ( * ) KK decays, the branching ratios are at the order from 10 −9 to 10 −7 . The two-body branching fraction B(B → D ( * ) R) can be extracted from the corresponding quasi-two-body decay modes in Tables II-V under the narrow width approximation relation By using the measured branching fractions B(φ(1020 [103] as an input, we can extract the branching ratios of two- , which agree with those from the two-body analyses based on the PQCD approach [23][24][25] within errors. The tiny differences between our predictions of the two-body decays and previous results of Refs. [23][24][25] are mainly due to the parametric origins and the power corrections related to the ratio r 2 As we know, a significant impact of nonfactorizable contribution is expected for a colour suppressed decay mode. Taking the decay channel B 0 s →D * 0 (φ(1020) →)K + K − as an example, one can see that the inclusion of the power correction r 2 D can suppress the branching ratio efficiently, especially for the contributions of nonfactorizable emission diagrams as shown in Table VII. Now, we come to discuss the contributions of P -wave resonances. It is obvious that the PQCD pre- are consistent well with the experimental data within errors. Taking B 0 s →D 0 K + K − decay as an example, we can calculate the total P -wave resonance contributions B(B 0 s →D 0 K + K − ) P-wave = 1.32 × 10 −5 , which is nearly equal to the sum of resonance contributions from φ(1020) and φ(1680). The main reason is that the interference between the two P -wave resonant states φ(1020) and φ(1680) is really small due to the rather narrow width of the former (Γ φ(1020) = 4.25MeV). Since the contribution from high mass resonance is one order of magnitude smaller, the P -wave resonance contribution is dominant by φ(1020). From the numerical results given in Table II and III, we obtain the relative ratio R(φ) between the branching ratio of B 0 s meson intoD * 0 and TABLE II: PQCD results for the branching ratios of the S, P and D wave resonance channels in the B 0 (s) →DK + K − decay together with experimental data [103]. The theoretical errors are attributed to the variation of the shape parameters ω B (s) (CD) in the wave function of B (s) (D) meson and the decay constant f B (s) , the Gegenbauer moments of two-kaon DAs, and the hard scale t and the QCD scale Λ QCD , respectively.

Decay Modes
Quasi-two-body  III: PQCD results for the branching ratios of the S, P and D wave resonance channels in the B 0 (s) → D * K + K − decay together with experimental data [103]. The theoretical errors are attributed to the variation of the shape parameters ω B (s) (CD * ) in the wave function of B (s) (D * ) meson and the decay constant f B (s) , the Gegenbauer moments of two-kaon DAs, and the hard scale t and the QCD scale Λ QCD , respectively.

Decay Modes
Quasi-two-body which is in agreement with the corresponding ratio of φ reported by the LHCb measurement [18], where the first uncertainty is statistical and the second is systematic. Utilizing the PQCD prediction in Table II and B(B 0 s →D 0K * 0 )) = (2.33 +1.18 −0.69 ) × 10 −4 taken from our previous work in Ref. [59] together with the two known branching ratios B(φ → K + K − ) = (49.2 ± 0.5)% V: PQCD results for the branching ratios of the S, P and D wave resonance channels in the B 0 (s) → D * K + K − decay. The theoretical errors are attributed to the variation of the shape parameters ω B (s) (C D * ) in the wave function of B (s) (D * ) meson and the decay constant f B (s) , the Gegenbauer moments of two-kaon DAs, and the hard scale t and the QCD scale Λ QCD , respectively.

Decay Modes
Quasi-two-body which is a bit larger than the data measured by LHCb [14], Recalling that the theoretical errors are relatively large, one still can count them as being consistent within one standard deviation. In contrast to the vector resonances, the identification of the scalar mesons is a long-standing puzzle. Scalar resonances are difficult to resolve because some of them have large decay widths, which cause a strong overlap between resonances and background. The prominent appearance of the f 0 (980) implies a dominant (ss) component in the semileptonic D s decays and decays of B (s) mesons. Ratios of decay rates of B and/or B s mesons into J/ψ plus f 0 (980) or f 0 (500) were proposed to allow for an extraction of the flavor mixing angle and to probe the tetraquark nature of those mesons within a certain model [125,126]. The phenomenological fits of the LHCb collaboration do neither allow for a contribution of the f 0 (980) in the B → J/ψππ [121] nor an f 0 (500) in B s → J/ψππ decays [120] by using the isobar model. The authors conclude that their data is incompatible with a model where f 0 (980) is formed from two quarks and two antiquarks (tetraquarks) at the eight standard deviation level. In addition, they extract an upper limit for the mixing angle of 17 • at 90% confidence level between the f 0 (980) and the f 0 (500) that would correspond to a substantial (ss) content in f 0 (980) [121]. However, a substantial f 0 (980) contribution is also found in the B-decays in a dispersive analysis of the same data that allows for a model-independent inclusion of the hadronic final state interactions in Ref. [127], which puts into question the conclusions of Ref. [121]. At this stage, the quark structure of scalar particles are still quite controversial. As a first approximation, the S-wave time-like form factor F S (m 2 KK ) used to parameterize the S-wave two-kaon distribution amplitude has been determined in Ref. [61]. Therefore, we take into account f 0 (980) and f 0 (1370) in thess density operator.
2 →]K + K − decays together with experimental data [103]. The theoretical errors are attributed to the variation of the shape parameters ω B (s) (C D * ) in the wave function of B (s) (D * ) meson and the decay constant f B (s) , the Gegenbauer moments of two-kaon DAs, and the hard scale t and the QCD scale Λ QCD , respectively. For the S-wave resonances f 0 (980) and f 0 (1370), firstly we define the ratio between the f 0 (980) → K + K − and f 0 (980) → π + π − ,
From the numerical results given in Table II, we can evaluate the relative branching ratios between two tensor modes, which can be tested by the forthcoming LHCb and Belle-II experiments. In Fig. 3, we show the ω-dependence of the differential decay rate dB(B 0 s →D 0 K + K − )/dω after the inclusion of the possible contributions from the resonant states f 0 (980) (the black solid curve), f 0 (1370) (the red solid curve), φ(1020) (the orange dotted curve), φ(1680) (the blue solid curve), f 2 (1270) (the cyan solid curve), f ′ 2 (1525) (the magenta solid curve), f 2 (1750) (the violet solid curve) and f 2 (1950) (the navy solid curve). To see clearly all the peaks of the various resonances, we draw them both in Fig. 3(a)  linear scale and Fig. 3(b) with a logarithmic one. For the considered decay mode B 0 s →D 0 K + K − , the dynamical limit on the value of invariant mass ω is . It is clear that an appreciable peak arising from the φ(1020) resonance, followed 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, while their broader widths compensate the integrated strength over the whole phase space. Therefore, the branching ratios of the four components are of a comparable size as predicted in our work. The contribution of the tensor f 2 (1270) is really small in B 0 s →D 0 K + K − decay mode, since it decays predominantly into ππ. Apart from above obvious signal peak, there are two visible structures at about 1750 MeV and 1950 MeV in Fig. 3(b). Obviously, the differential branching ratios of these decays exhibit peaks at the pole mass of the resonant states. Thus, the main portion of the branching ratio lies in the region around the resonance as expected. For B 0 s →D 0 φ →D 0 K + K − decay, the central values of the branching ratio B are 5.82 × 10 −6 and 8.29 × 10 −6 when the integration over ω is limited in the range of respectively, which amount to 48.1% and 68.5% of the total branching ratio B = 1.21 × 10 −5 as listed in Table II.
2 →]K + K − decays are vector-vector (vector-tensor) modes and can proceed through different polarization amplitudes. We display the polarization fractions associated with the available data in Table VI, which have the same origin of theoretical uncertainties as the branching ratios. It is easy to see that the fraction of the longitudinal polarization can be generally reduced to about ∼ 50%, while the parallel and perpendicular ones are roughly equal. The results are quite different from the expectation in the factorization assumption that the longitudinal polarization should dominate based on the quark helicity analysis [132,133]. For B 0 s →D * 0 (φ(1020) →)K + K − decay, although the central value of f 0 ≈ 55% is a little smaller than the measured one f exp 0 = 73%, but they are consistent with each other due to the still large theoretical and experimental uncertainties.
It is worthwhile to stress that the polarization fractions of the colour-suppressed decays B →D * (f 2 → )KK and B → D * (f 2 →)KK are quite different from each other. For B →D * (f 2 →)KK (b →c transition) decays, the percentage of the transverse polarizations (f + f ⊥ ) can be as large as 80%, while for B → D * (f 2 →)KK (b →ū transition) decays, the percentage of the transverse polarizations are only at the range of 10% ∼ 20%. This situation can be understood as the following reasons [25]: TheD * meson in B →D * (f 2 →)KK decays is longitudinal polarized since thec quark and the u quark in theD * meson produced through the (V −A) current are right-handed and left-handed, respectively. Because thec quack is massive, the helicity of thec quark can flip easily from right handed to left handed. Thus, theD * meson can be transversely polarized with the polarization λ = −1. The recoiled tensor meson can also be transversely polarized with the polarization λ = −1 due to the contribution of orbital angular momentum. Then the transverse polarized contribution in B →D * (f 2 →)KK decays can be sizable. While for the D * meson, theū quark in the D * meson is right handed and the c quark can flip from left handed to right handed, which makes the D * meson transversely polarized with the polarization λ = +1. In order to have the same polarization with the D * meson, the recoiled transversely polarized tensor meson needs contributions from both orbital angular momentum and spin, so the situation is symmetric. But the wave function of the tensor meson is asymmetric. Therefore, the transversely polarized contribution in B → D * (f 2 →)KK is suppressed on account of Bose statistics. More precise measurements of such decay channels are expected to help us to test and improve our theoretical calculations. Each partial wave contribution can be parameterized into the corresponding time-like form factor, which contains the final state interactions between the kaons in the resonant regions. The Flatté model for the f 0 (980) resonance and the Breit-Wigner formulas for other resonances have been adopted to parameterize the time-like form factors F S,P,D (ω 2 ) involved in the dikaon DAs. Using the determined parameters of the two-kaon DAs in our previous work, we have predicted the branching ratios of B (s) → [D ( * ) ,D ( * ) ](R → )KK decay channels. It has been shown that our predictions of the branching ratios for most of the considered decays are in good agreement with the existing data within the errors. The branching ratios of the twobody B (s) → [D ( * ) ,D ( * ) ]R can be extracted from the corresponding quasi-two-body modes by employing the narrow width approximation. Moreover, we calculated the polarization fractions of the vector-vector and vector-tensor decay modes in detail. For most of the considered channels, the transverse polarization has been found to be of the similar size as the longitudinal one, which is quite different from the general expectation in the factorization assumption. More precise data from the LHCb and the future Belle II will test our predictions.

Appendix A: Decay amplitudes
For B (s) →D ( * ) (R →)KK decays (R denotes the various partial wave resonances ), the expressions of the individual amplitudes F LL e(a) and M LL e(a) can be straightforwardly obtained by evaluating the Feyman diagrams in Fig. 1. Performing the standard PQCD calculations, one gets the following expressions of the relevant amplitudes: with the ratio r D = m D ( * ) /m B (s) and the color factor The formulas for the longitudinal component amplitudes A L are as follows: The transverse polarization amplitudes A N,T are of the following form: M LL eφ,T = 16 M LL aφ,N = 16 M LL aφ,T = 16 For the decays involved D ( * ) mesons, similarly, the corresponding decay amplitudes can be obtained by evaluating the Feyman diagrams in Fig. 2.
The hard scales t i appeared in the above equations are chosen as the maximum of the virtuality of the internal momentum transition in the hard amplitudes. For B (s) →D ( * ) (R →)KK decays, we have