Four-body decays $B_{(s)} \rightarrow (K\pi)_{S/P} (K\pi)_{S/P}$ in the perturbative QCD approach

We investigate the four-body decays $B_{(s)} \rightarrow (K\pi)_{S/P} (K\pi)_{S/P}$ in the $K\pi$-pair invariant mass region around the $K^*(892)$ resonance in the perturbative QCD (PQCD) approach, where $(K\pi)_{S/P}$ denotes a $S$- or $P$-wave $K\pi$ configuration. The contributions from the $P$-wave resonance $K^*(892)$ and from possible $S$-wave components are included into the parametrization of the $K\pi$ two-meson distribution amplitudes, which are the nonperturbative inputs in this formalism. We calculate the branching ratio for each resonance and observe sizable $S$-wave contributions to the $B_s$ modes, well consistent with recently available LHCb data. We also deduce the polarization fractions for the $K^*\bar{K}^*$ final states, and compare our predictions to those in previous PQCD analyses of the corresponding two-body $B\to K^*\bar{K}^*$ decays. The polarization puzzle associated with the two $U$-spin related channels $B^0\rightarrow K^{*0}\bar{K}^{*0}$ and $B^0_s\rightarrow K^{*0}\bar{K}^{*0}$ still exists. In addition to direct $CP$ asymmetries, triple-product asymmetries and $S$-wave induced direct $CP$ asymmetries originating from the interference among various helicity amplitudes, are presented. It is shown that true triple-product asymmetries are rather small, while direct $CP$ asymmetries and $S$-wave-induced direct $CP$ asymmetries are significant in some decays. Our results will be subject to stringent tests with precise data from $B$ factories in the near future.


I. INTRODUCTION
The charge conjugation (C), parity (P ), and time reversal (T ) are three fundamental discrete transformations in nature. The combination of the C and P transformations is particularly interesting, because CP violation is a key requirement for generating the matter-antimatter asymmetry in the universe [1,2]. The presence of CP violation, attributed to an irreducible phase in the Cabibbo-Kobayashi-Maskawa (CKM) quark-mixing matrix [3,4] in the Standard Model (SM), has been well established experimentally in the K and B meson systems [5]. CP violation has been observed not only in two-body B meson decays, but also in multi-body channels. In addition to direct and mixing-induced CP asymmetries, triple-product asymmetries (TPAs) [6][7][8] can be defined in multi-body B meson decays, which may reveal potential signals of CP violation. The construction of a scalar triple product requires a final state with at least four particles, which endow three independent momenta in the rest frame of a B meson. That is, four-body B meson decays are rich in CP violation phenomena in the quark sector, and the extra asymmetry observables provide more opportunities to search for new physics beyond the SM.
The B → V 1 (→ P 1 P ′ 1 )V 2 (→ P 2 P ′ 2 ) decay, where four pseudoscalar mesons are produced through two intermediate vector resonances, gives rise to three helicity amplitudes commonly referred to as the P -wave amplitudes. There may exist a background from the decay of a scalar resonance S → P P ′ or from a scalar nonresonant P P ′ production [9,10], such that three more amplitudes contribute: two single S-wave amplitudes (the scalar and vector combination) and one scalar-scalar double S-wave amplitude. The interference among the above helicity amplitudes generate TPAs, which can be extracted from an angular analysis for the decay. A scalar triple product takes a generic form , v i being the spin or momentum of a final-state particle. This triple product is odd under the time reversal transformation, and also under the CP transformation due to the CP T invariance. A TPA A T represents an asymmetry between the decay rates Γ with positive and negative values of T P , A T ≡ Γ(T P > 0) − Γ(T P < 0) Γ(T P > 0) + Γ(T P < 0) .
In the presence of the scalar background, S-wave induced direct CP asymmetries, besides TPAs, are accessible in untagged decays [11]. Phenomenological investigations on these asymmetry observables have been conducted intensively in the literature [11][12][13][14][15][16]. B meson decays into various four-body charmless hadronic final states in certain two-body invariant mass regions have been observed by the Belle [17][18][19], BaBar [20][21][22], and LHCb [23][24][25][26][27][28][29][30] Collaborations. Branching ratios for many partial waves have been measured for the first time or with greatly improved accuracy. The B 0 s → (K + π − )(K − π + ) mode was first seen by the LHCb [23], whose signals mainly come from K * K * channels with some S-wave contributions, as indicated by Kπ mass distributions. The involved amplitudes were determined later for the K + π − and K − π + invariant masses around the K * (892) resonance in the window of ±150 MeV [26]. The polarization fractions, TPAs, and S-wave fractions were also extracted through a combined angular and mass analysis. Subsequently, a flavortagged decay-time-dependent amplitude study of this mode in a K ± π ∓ mass window extending from 750 to 1600 MeV was presented [28], where not only the K * andK * resonances, but also several scalar, vector and tensor components in the Kπ system were examined. More recently, the LHCb performed a combined amplitude analysis for the B 0 s → (K + π − )(K − π + ) and B 0 → (K + π − )(K − π + ) decays, which are related by the U -spin symmetry, and found different polarization patterns: the latter is dominated by the longitudinal polarization, while the low longitudinal polarization contribution to the former deviates from theoretical expectations in the SM [31][32][33][34][35][36][37].
Recently, the localized CP violation and branching fraction of the four-body decayB 0 → K − π + π + π − were calculated by employing a quasi-two-body QCDF approach in Refs. [72,73]. Similar to the handling of three-body B meson decays in the PQCD approach, four-body processes are assumed to proceed dominantly with two intermediate resonances, which then strongly decay into two light meson pairs. The PQCD factorization formalism for four-body B meson decays is thus simplified to the one for two-body decays [74,75]. A decay amplitude is then written as where Φ B is the B meson DA, and the two-meson DAs Φ P1P ′ 1 and Φ P2P ′ 2 absorb the nonperturbative dynamics involved in the meson pairs P 1 P ′ 1 and P 2 P ′ 2 , respectively. The corresponding hard kernel H is evaluated at the quark level in perturbation theory. The symbol ⊗ stands for the convolution of all the perturbative and nonperturbative factors in parton momenta.
In the present work we focus on the four-body decays B (s) → (Kπ) S/P (Kπ) S/P with (Kπ) S/P denoting a S-or P -wave Kπ configuration, and derive their decay amplitudes in the PQCD approach. Six quasi-two-body channels are computed, which arise from different combinations of the Kπ pairs in S and P waves. The Kπ two-meson DAs have been constrained to some extent in previous PQCD studies of three-body B meson decays [76][77][78]. With these universal nonperturbative inputs, we make quantitative predictions for physical observables in the above modes, including branching ratios, polarization fractions, S-wave fractions, direct CP asymmetries, and TPAs. Note that a nonzero TPA can be generated by either final-state interactions or CP violation. To justify a "true" CP violation signal, we compare a TPA with that in the corresponding CP conjugate channel. In accordance with the LHCb measurements [26,30], the invariant masses of the Kπ pairs are restricted to the region around the K * (892) resonance in the window of 150 MeV. The formalism presented here, ready to be extended to other four-body hadronic B meson decays, will have many potential applications.
The paper is organized as below. In Sec. II we define the kinematic variables for four-body B meson decays and the two-meson DAs for the S-and P -wave Kπ configurations. The B (s) → (Kπ) S/P (Kπ) S/P helicity amplitudes are formulated in Sec. III, in terms of which the polarization fractions and direct CP asymmetries are expressed. The various TPAs and S-wave-induced direct CP asymmetries to be investigated are specified in Sec. IV. Numerical results are presented and discussed in Sec. V, which is followed by the conclusion in Sec. VI. The explicit PQCD factorization formulas of all the helicity amplitudes are collected in Appendix A. Consider the B meson decay into two intermediate states M 1 and M 2 , which further strongly decay into the pseudoscalar pairs P 1 P ′ 1 and P 2 P ′ 2 , respectively, with the meson momenta p B = p + q, p = p 1 + p ′ 1 , and q = p 2 + p ′ 2 . We work in the rest frame of the B meson, and choose the momentum p B = (M/ √ 2)(1, 1, 0 T ) in the light-cone coordinates, M being the B meson mass. The momenta of the two meson pairs can be parametrized as where the factors with the mass ratios η 1,2 = ω 2 1,2 /M 2 , are related to the invariant masses of the meson pairs via p 2 = ω 2 1 , q 2 = ω 2 2 , and p B = p + q. For the P -wave pairs, the corresponding longitudinal polarization vectors are defined as which satisfy the normalization ǫ 2 p = ǫ 2 q = −1 and the orthogonality ǫ p · p = ǫ q · q = 0. We derive the momentum of each final-state meson, from the relations p = p 1 + p ′ 1 and q = p 2 + p ′ 2 , and the on-shell conditions p (4) and (7) give the meson momentum fractions where the variables ζ i bear the meaning of the meson momentum fractions up to corrections from the meson masses. The transverse momenta p T and q T are written as with the factors Alternatively, one can define the polar angles θ i of the mesons P i in the P i P ′ i rest frames and the azimuthal angle φ between the decay planes of M 1 and M 2 as illustrated in Fig. 1 to describe the sequential decay. The transformation connecting the B meson rest frame and the meson pair rest frame leads to the relations between the meson momentum fractions ζ i and the helicity angles θ i , with the bounds Note that Eq. (11) reduces to the conventional expression in [79,80] when the two mesons in a pair have equal masses, and that it maintains the decoupling between S-and P -wave contributions to a total decay rate in the case with arbitrary final-state meson masses: We emphasize that the parametrization with the exact dependence on the final-state meson masses in Eq. (7) is crucial for establishing Eq. (11), and then Eq. (13). The Feynman diagrams contributing to the hard kernels for the B → M 1 M 2 → P 1 P ′ 1 P 2 P ′ 2 decays in the PQCD approach are displayed in Fig. 2, each of which contains a single virtual gluon exchange at leading order in the strong coupling α s . The diagrams in the first (second) row correspond to the emission (annihilation) type, which are further classified into the factorizable ones with hard gluons attaching to the quarks in the same mesons, and the nonfactorizable ones with hard gluons attaching to the quarks in different mesons. To evaluate the hard kernels, we parametrize the three valence quark momenta labelled by k B , k 1 , and k 2 in Fig. 2(a) as with the parton momentum fractions x i and the parton transverse momenta k iT , i = B, 1, 2. Since k 1 (k 2 ) is aligned with the meson pair in the plus (minus) direction, its small minus (plus) component has been neglected. We also drop where the symbol • denotes a weak vertex, and the cyan oval represents either a scalar or vector resonance.
the plus component k + B , because it does not appear in the hard kernels for dominant factorizable contributions. The next-to-leading-order corrections to the hard kernels, such as those from quark loops [81], will be taken into account in the future.
The light-cone hadronic matrix element for a B meson is decomposed as [82] with b being the impact parameter conjugate to the parton transverse momentum k T , and N c being the number of colors. A model-independent determination of the B meson DA from an Euclidean lattice was attempted very recently [83]. Various models of φ B are available in the literature [84,85], among which we adopt the conventional one [82,86], with the shape parameter ω b = 0.40 GeV for B u,d mesons and ω b = 0.48 GeV for a B s meson [87]. The normalization constant N B is related to the B meson decay constant f B via the normalization Note that the above model φ B (x, b = 0) corresponds to the inverse moment λ B = 241 (λ Bs = 287) MeV for the B (B s ) meson with the mass M B = 5.28 (M Bs = 5.37) GeV. This λ B value is lower than 460 ± 110 MeV evaluated at 1 GeV scale in QCD sum rules [88], but consistent with the experimental bound λ B > 238 MeV [89]. The two-and three-parton B meson DAs have been decomposed according to definite twist and conformal spin assignments up to twist 6 in Ref. [90]. Here only the B meson DA associated with the leading Lorentz structure in Eq. (15) is considered, while other power-suppressed pieces are negligible within the accuracy of the current work [86,91]. It has been shown explicitly in the PQCD approach that the next-to-leading-power B meson DA contributes about 17% of the B → π transition form factors at large recoil [92], consistent with the naive estimate for the power suppression effectΛ/m b ≡ (M B − m b )/m b ∼ 10%, m b being the b quark mass. We will include the power-suppressed B meson DAs, together with other subleading contributions [93,94], for a complete analysis of four-body B meson decays in the future, when aiming at precision better than 10%.
As stated before, we focus on the leading-power regions of phase space, where the four-body decays B (s) → (Kπ) S/P (Kπ) S/P proceed mainly via quasi-two-body processes. Besides the B meson DA, the two-meson DAs, which absorb strong interaction related to the production of the Kπ meson pairs, are the necessary inputs to PQCD calculations. In what follows the subscripts S and P label the corresponding partial waves. The light-cone matrix element for a S-wave Kπ pair in the plus direction is decomposed as [76,78] with the dimensionless vectors n = (1, 0, 0 T ) and v = (0, 1, 0 T ), and the Kπ invariant mass ω. The explicit forms of the twist-2 DA φ S and twist-3 DAs φ s S and φ t S are given by [76]: with the ratio µ S = ω/(m s − m q ), where m s (1 GeV) = 119 MeV [95,96] is the running strange quark mass and the light quark masses m q , q = u, d, are set to zero. We take the Gegenbauer moments B 1 = −0.57 ± 0.13 and B 3 = −0.42 ± 0.22 at the 1 GeV scale from scenario I in the QCD sum rule analysis [95,96], which have been widely employed in the PQCD studies of the B → K * 0 (1430) transitions [97][98][99]. Another possible choice of the Gegenbauer moments as those for the κ(800) meson DA does not work well for accommodating experimental data. Because available data are not yet sufficiently precise for fixing more Gegenbauer moments, the asymptotic forms have been assumed for the two twist-3 DAs.
For the scalar form factor F s (ω 2 ), we follow the LASS line shape [100], which consists of the K * 0 (1430) resonance as well as an effective-range nonresonant component, with the K * 0 (1430) mass (width) m 0 = 1450 ± 80 MeV (Γ 0 = 400 ± 230 MeV) [101]. The kaon three-momentum k(ω) is written, in the Kπ center-of-mass frame, as m K(π) being the kaon (pion) mass. For the scattering length a and effective-range parameters b, we take the values a = 3.2 ± 1.8 GeV −1 and b = 0.9 ± 1.1 GeV −1 [101]. The light-cone matrix elements for longitudinal and transverse P -wave Kπ pairs are organized, in analogy with the dipion case [102], into respectively. Note that the factor 2ζ − 1 in Φ L P corresponds exactly to the structure for a time-like vector form factor, with the kinematic variables in Eqs. (4), (6) and (7) being inserted. This is why the first Lorentz structure for the longitudinal Kπ pair has been set to the polarization vector ǫ p , different from the parametrization in [77]. The factor ζ(1 − ζ) + α for the transverse Kπ pair comes from Eq. (9).
The twist-2 DAs φ P and φ T P , and the other twist-3 DAs appearing in Eq. (22) are expanded in terms of the Gegenbauer polynomials [103], with the Gegenbauer moments [103] That is, we do not distinguish the Gegenbauer moments for the longitudinal and transverse polarizations. The time-like vector form factor F Kπ is parametrized as the relativistic Breit-Wigner (RBW) function [77,78] with the mass-dependent width where m K * (Γ K * ) is the K * (892) mass (width), and r = 3.0 GeV −1 [26] is the interaction radius. Here only the K * (892) resonance is included, and other higher resonances, such as K * (1410) and K * (1680) with pole masses much above the considered mass window, are neglected. We point out that a phase difference causing the interference between the S-and P -wave form factors is expected, which may be a complicated function of the Kπ invariant mass. Concentrating on the region near the K * resonance, we perform the Taylor expansion of this complicated function around m K * , and keep only the leading nontrivial term proportional to m 2 K * − ω 2 . The first term in the expansion, giving rise to an overall constant phase, does not contribute to the interference effect. The vector form factor is thus replaced by 1 where the tunable parameter β renders the phase dimensionless. As shown in Sec. V, the choice β = (7.5 ± 2.5) GeV −2 improves the agreement between theoretical results and data. For another form factor F ⊥ Kπ , we assume the with the tensor (vector) decay constant of the K * (892) meson f T K * = 0.185 GeV (f K * = 0.217 GeV) [105] derived in the narrow-width limit in QCD sum rules. Note that the tensor decay constant is renormalization-scale and renormalization-scheme dependent, whose value is taken at the typical scale 1 GeV, and in the M S scheme. We assume the naive relation between the Kπ time-like form factors F Kπ and F ⊥ Kπ , due to lack of rigorous theoretical and experimental information. For example, the establishments of the QCD sum rules for the form factors depend on their parametrization [106], and the determinations of the decay constants are altered by finite-width effects from a K * meson [106].

III. HELICITY AMPLITUDES
The differential rate for the decays B (s) → (Kπ) S/P (Kπ) S/P in the B (s) meson rest frame is written as where dΩ with Ω ≡ {θ 1 , θ 2 , φ, ω 1 , ω 2 } stands for the five-dimensional measure spanned by the three helicity angles and the two invariant masses, and is the momentum of one of the Kπ pairs in the B (s) meson rest frame. The four-body phase space has been derived in the analyses of the K → ππlν decay [107], the semileptonicB → D(D * )πlν decays [108], semileptonic baryonic decays [109,110], and four-body baryonic decays [111]. One can confirm that Eq. (34) is equivalent to those in Refs. [109,111] by appropriate variable changes. Replacing the helicity angles θ i by the meson momentum fractions ζ i via Eq. (11), we turn Eq. (34) into The B (s) → (Kπ) S/P (Kπ) S/P decays involve six helicity amplitudes A h with h = V V (three), V S, SV , and SS. The first three, commonly referred to as the P -wave amplitudes, correspond to the final states, where both Kπ pairs are from intermediate vector mesons. A P -wave decay amplitude can be decomposed into three components in the transversity basis: A 0 , for which the polarizations of the vector mesons are longitudinal with respect to their momenta, and A (A ⊥ ), for which the polarizations are transverse to the momenta and parallel (perpendicular) to each other. As the S-wave Kπ pair arises from M 1 (M 2 ) in Fig. 2(a), the resultant single S-wave amplitude is denoted as A SV (A V S ). The double S-wave amplitude A SS is associated with the final state, where both Kπ pairs are produced in the S wave. Specifically, these helicity amplitudes represent the following B 0 where the notation (Kπ) 0 labels the S-wave Kπ configuration to be modelled by the LASS line shape. It is more convenient to introduce the linear combinations for the discussion of asymmetry observables below. Note that the amplitude A SV (A V S ) in our convention corresponds to A V S (A SV ) in Ref. [26], so that the above definitions of A S ± are the same as theirs. Note that the amplitudes A ⊥ and A S + are CP -odd, and A 0 , A , A S − and A SS are CP -even for the neutral B 0 (s) modes. Including the ζ i and azimuth-angle dependencies, we have the full decay amplitude in Eq. (36) as a coherent sum of the P -, S-, and double S-wave components, To calculate the helicity amplitudes, we start from the effective Hamiltonian [112] with the quark q = d, s, the Fermi coupling constant G F , the CKM matrix elements V ub , V uq , · · · , and the Wilson coefficients C i (µ) at the renormalization scale µ. The effective local four-quark operators O i (µ) are given by [112] where α, β are color indices and e q ′ is the electric charges of the quark q ′ in units of |e|. The sum over q ′ runs over the quark fields active at the b quark mass scale, i.e., q ′ = u, d, s, c, b.
As stated in the Introduction, a helicity amplitude is expressed as a convolution of a hard kernel with relevant meson DAs in the PQCD approach. Parton transverse momenta in the k T factorization theorem, on which this approach is based, regularize end-point singularities from the kinematic regions with small momentum fractions x. The resultant double logarithms α s ln 2 k T are organized into the Sukakov factor e −S through the k T resummation technique. The double logarithms α s ln 2 x appearing in a hard kernel, being important in the end-point regions, are summed to all orders into the threshold resummation factor S t . The typical PQCD factorization formula then reads where b i are the impact parameters conjugate to k iT in Eq. (14), t represents the largest energy scale in the hard kernel H, and "Tr" denotes the trace over all Dirac structures and color indices. Both the resummaton factors e −S and S t , whose explicit expressions are provided in Appendix A, improve the perturbative evaluation of a helicity amplitude. The prefactors compensated by the denominators √ 1 + 4α i in Eq. (39), normalize the helicity amplitudes A h properly.
The helicity amplitudes for the pure penguin type channel B 0 s → K * 0 (→ K − π + )K * 0 (→ K + π − ) are written as where and (S − P ) (S + P ) operators, respectively. The explicit expressions of all the above functions can be found in Appendix A. Note that the function F SP e from the operators O 5−8 vanishes, when a vector resonance is emitted from the weak vertex, because neither the scalar nor pseudoscalar operator is responsible for vector resonance production. Equation (44) also holds for the B 0 → K * 0 (→ K − π + )K * 0 (→ K + π − ) decay with the replacement V * tb V ts → V * tb V td . Additional tree contributions appear in the B 0 (s) → K * + (→ K 0 π + )K * − (→K 0 π − ) and B + → K * + (→ K 0 π + )K * 0 (→ K + π − ) decays, whose helicity amplitudes are given by The isospin ratio Γ(K * + → K + π 0 )/Γ(K * + → K 0 π + ) = 1/2 leads to the relations We obtain the branching ratios from Eq. (36), where the invariant masses ω 1,2 are integrated over the selected Kπ mass window, and the coefficients are the results of the integrations over ζ 1 , ζ 2 , φ. In terms of Eq. (49), we compute the CP -averaged branching ratio and the direct CP asymmetry in each component, respectively, whereB h is the branching ratio of the corresponding CP -conjugate channel. The sum of the six components yields the total branching ratio and the overall direct CP asymmetry, respectively.
To characterize the relative importance of the S-wave contributions, we define the S-wave fractions for the different components as and the total S-wave fraction as The polarization fractions from the P -wave amplitudes are derived according to with B P = B 0 + B + B ⊥ being the total P -wave branching ratio.

IV. TRIPLE PRODUCT ASYMMETRIES AND S-WAVE-INDUCED DIRECT CP ASYMMETRIES
In this work we focus on the TPAs associated with A ⊥ for the B (s) → (Kπ) S/P (Kπ) S/P decays, which are defined in terms of the partially integrated differential decay rates as [12,26] with the denominator The above TPAs contain the integrands Im(A ⊥ A * h ) = |A ⊥ ||A * h | sin(∆φ+ ∆δ), where ∆φ and ∆δ represent the weak and strong phase differences between the amplitudes A ⊥ and A h , respectively. Note that a strong phase difference yields a TPA even in the absence of weak phases, so a nonzero TPA does not necessarily signal CP violation. To identify a true CP violation signal, one has to compare the TPAs in B andB meson decays. The helicity amplitudē A h , whose weak phases flip sign relative to those in A h , and the TPAsĀ i T for a CP -conjugated process can be derived similarly. We then have the true and fake asymmetries relative to the CP -averaged decay rate [12] respectively, withΓ being the decay rate of the CP -conjugate process. Note that A i T (true), nonvanishing only in the presence of the weak phase difference, provides an alternative measure of CP violation. Furthermore, compared with direct CP asymmetries, A i T (true) does not suffer the suppression from the strong phase difference, and reaches a maximum when the strong phase difference vanishes. On the contrary, A i T (f ake) can be nonzero even as the weak phase difference vanishes. Such a quantity is sometimes referred to as a fake asymmetry, which reflects the effect of strong phases [12], instead of CP violation. Since the helicity amplitudes may have different strong phases, it is likely that fake TPAs do not diminish for all channels. In the modes with the neutral intermediate states K * 0K * 0 , each helicity amplitude involves the same single weak phase in the SM, such that true TPAs vanish without the weak phase difference. Therefore, true TPAs in these modes, if observed, are probably signals of new physics.
The interference of the amplitude A S + with the others generates more asymmetries [26]: Similarly, combining the quantities for the corresponding CP -conjugate process, we define the S-wave-induced direct CP asymmetries [26] A Note that the above quantities, as the sum of the asymmetries in the B andB meson decays, can be measured using untagged samples, in which CP -conjugate processes need not be distinguished [15].

V. NUMERICAL RESULTS
In this section we calculate the physical observables for the B (s) → (Kπ) S/P (Kπ) S/P decays raised in the previous sections, including the CP averaged branching ratios, S-wave fractions, polarization fractions, TPAs, and direct CP asymmetries. The inputs for the numerical analysis are summarized below. The meson masses are taken as (in units of GeV) [5] M Bs = 5.37, M B = 5.28, m K = 0.494, m π = 0.14.
For the CKM matrix elements, we adopt the Wolfenstein parametrization with the parameters [5] The lifetime (in units of ps) and decay constant (in units of GeV) of the B (s) meson are set to the values [5,31], Other parameters appearing in the Kπ two-meson DAs have been specified before. We first examine the B → Kπ transition form factors, to which the factorizable emission contributions from Figs. 2(a),(b) are related. These form factors serve as the nonperturbative inputs to QCDF evaluations of the factorizable emission amplitudes [47]. There are in general seven form factors for the P -wave B → Kπ transition and three for the S-wave one, whose definitions can be found in Refs. [106,113,114]. Here we focus only on the form factor F ⊥ (ω 2 1 , q 2 ) defined via the matrix element Kπ|sγ µ b|B with the momentum transfer q being given in Eq. (4). A complete and systematic PQCD investigation of all the B → Kπ form factors, including their ω 2 1 and q 2 dependencies, will be postponed to a future work. The corresponding PQCD factorization formula for F ⊥ (ω 2 1 , q 2 ) is also presented in Appendix A, which yields F ⊥ (m 2 K * (892) , 0) = 88 (64) for the B (s) → Kπ transition at maximal recoil. Results for the B → Kπ transition from light-cone sum rules in the narrow-width limit, ranging between 62 and 93 [106,[115][116][117], agree well with ours. Our B s → Kπ transition form factor is also consistent with 72 ± 7 derived in Ref [116].
A. CP averaged four-body branching ratios and S-wave fractions The integral in Eq. (49) over ω 1,2 in the selected mass range gives the CP -averaged branching ratios of various components and their sum for five penguin-dominated channels, which are given in Table I. The two B s modes, induced by the b → s transition, have relatively large branching ratios of O(10 −6 ). The remaining three modes, mediated by the b → d transition, are generally an order of magnitude smaller. In particular, the B 0 → (K 0 π + )(K 0 π − ) decay, proceeding only through the power-suppressed weak annihilation, has an even lower rate. We have considered three sources of theoretical uncertainties: the shape parameters in the B (s) meson DAs with 10% variation, the Gegenbauer moments in the S-and P -wave Kπ two-meson DAs, and the hard scales t defined in Appendix A, that vary from 0.75t to 1.25t. Their effects are listed in Table I in order, among which the third one is more dominant. Higher-order contributions to four-body B meson decays can be included to reduce the sensitivity to the variation of the hard scales. The Gegenbauer moments also need to be constrained for further improving the precision of theoretical predictions.
It is straightforward to derive from the branching ratios in Table I and Eq. (53) the S-wave fractions in each mode as presented in Table II. According to Eq. (38), the amplitudes A S − and A S + , as the linear superposition of A SV and A V S , exhibit the interference between the S-and P -wave configurations. Therefore, the phenomenological parameter β introduced in Eq. (33) affects the relative strength of f S + and f S − . Figure 3 illustrates the influence of β on the two single S-wave fractions, where the dashed red and solid blue curves describe f S + and f S − , respectively. It is found that the S-wave fractions are more sensitive to β in the region of [-5,20] GeV −2 . The value of β can be bounded in the range β = (7.5 ± 2.5) GeV −2 , in which f S + and f S − have reasonable relative magnitudes to match the data better. The component f SS , which is proportional to the modulus squared of its corresponding amplitude, is independent of β. The total S-wave contribution to f S-wave is also independent of the phase parameter, since the interference terms cancel in the sum f S + + f S − in Eq. (54).
The obtained three S-wave fractions for the B 0 s → (K + π − )(K − π + ) decay are compatible with the two measurements in Refs. [26,30] within errors. The central values of f SS and f S + for the B 0 → (K + π − )(K − π + ) decay are obviously larger, while f S − is a bit smaller than the data [30]. Nevertheless, the total S-wave contribution, accounting for about 42.4% of the total decay rate, is consistent with the data [30]. More precise measurements and more II: S-wave fractions in the B (s) → (Kπ) S/P (Kπ) S/P decays within the Kπ invariant mass window of 150 MeV around the K * (892) resonance. The data corresponding to the same invariant mass window are taken from [26,30], where the first uncertainty is statistical and the second is systematic.
complete theoretical analyses at the subleading level will help clarify the discrepancy. We point out that the S-wave contributions could be modified by adding higher Gegenbauer terms to the twist-3 DAs for the S-wave Kπ pair, and better consistency with the data is expected. However, such fine tuning will not be attempted in the present work. The S-wave fractions of the other channels in Table II exhibit similar tendency. The predicted total S-wave fractions, ranging from 29.5% to 51.4%, can be tested at future experiments. It is obvious that the S-wave contributions to the B (s) → (Kπ) S/P (Kπ) S/P decays in the considered mass region are significant.

−21
· · · · · · Our predicted branching ratio for the pure penguin mode B 0 s → K * 0K * 0 , being higher than those from the PQCD and QCDF approaches, agrees with the data. As to the branching ratio of another pure penguin mode B 0 → K * 0K * 0 , our prediction is slightly lower than the PQCD and QCDF ones, which are all systematically below the data. It has been found [118] that the factorizable emission amplitude F e gives the dominant contribution to the B 0 → K * 0K * 0 decay, and the annihilation amplitudes cancel strongly between F a /M a and F ′ a /M ′ a in Eq. (44). On the contrary, such cancellation does not occur in the charged mode B + → K * +K * 0 owing to vanishing F ′ a /M ′ a . Moreover, the charged mode receives a color-allowed tree contribution from the annihilation diagrams, which is induced by the (V −A)(V −A) operators with the large Wilson coefficient C 2 + C 1 /3, but subject to the helicity suppression compared with that from the emission diagrams. The interference between the annihilation and emission amplitudes turns out to increase the branching ratio of the charged mode, such that our branching ratio B(B + → K * +K * 0 ) is about twice of the neutral one B(B 0 → K * 0K * 0 ). It is thus greater than those from the other approaches and more consistent with the data. The pattern B(B + → K * +K * 0 ) > B(B 0 → K * 0K * 0 ) has been also observed in the previous PQCD investigations [31,118].
The derived B 0 s → K * + K * − branching ratio is close to the B 0 s → K * 0K * 0 one, because the additional tree contribution to the former is minor. This feature is common among the predictions from all the approaches. The predicted branching ratio for the pure annihilation decay B 0 → K * + K * − spans a wide range, (0.064 − 1.43) × 10 −6 , with large theoretical uncertainties, as indicated in Table III. Our result is basically the same as the previous PQCD calculation [31], and about twice the QCDF one [119,120]. The value from the FAT approach [37] is far too large, but still lower than the upper experimental bound [5]. More precise measurements of this channel are crucial for testing the theoretical expectations.
The heavy quark limit implies that charmless B → V V decays mainly produce longitudinally polarized final states due to the vector-axial structure of weak currents and the quark helicity conservation in high-energy QCD interactions [35,121]. However, data for penguin-dominated B meson decays, contradicting to this speculation in general, hint that the dynamics of penguin transitions are more complicated. This so-called polarization puzzle in penguin-dominated B → V V modes has posed a challenge to the development of theoretical frameworks for B meson decays.
As shown in Table III, the longitudinal polarization fraction f 0 for the B 0(+) → K * 0(+)K * 0 decays from the PQCD approach (including the present work) are higher than the transverse one f T = f +f ⊥ . Both the QCDF [119,120] and SCET [36] yield the similar pattern f 0 ∼ f T , albeit with large uncertainties. The world average data in Table III [21,30,122,123] Table III that various predictions for f 0 (B 0 s → K * 0K * 0 ) lie in the range 0.343 to 0.636, considerably larger than the data. The QCDF evaluation gives a low longitudinal polarization fraction f 0 (B 0 s → K * 0K * 0 ) = (27.7 +8.2+9.5 −6.7−18.9 )% [124] by including weak annihilation corrections with the best-fit endpoint parameters. This abnormal longitudinal polarization fraction seems to be reconciled, but an obvious tension between the data and the prediction for the branching ratio B(B s → φK * ) is invoked accordingly. The PQCD approach yields a sizable transverse polarization, which is enhanced by the weak annihilation from the operator O 6 and by nonfactorizable contributions [125]. However, these effects are not able to fully account for the above polarization anomaly. The small f 0 = 0.38 in Ref. [31] is ascribed to the inclusion of the higher-power terms proportional to r 2 , with r being the mass ratio between the vector and B mesons. Our predictions for the longitudinal polarization fractions agree with the QCDF ones [34,35], and a bit higher than the previous PQCD results. It is worth mentioning that we have employed the Gegenbauer moments for the transversely polarized Kπ DAs the same as for the longitudinal polarized ones (see Eq. (30)) in this analysis. It should not be difficult to accommodate both the branching ratio and polarization data for the B s modes simultaneously by tuning the Gegenbauer moments for the former.
It is noted that the longitudinal polarization contributions to the B s → K * 0K * 0 and B 0 → K * 0K * 0 decays can be related to each other in the U -spin symmetry. Nevertheless, the data reveal quite different longitudinal polarization contributions in these two modes. This discrepancy indicates that the U -spin relations are not well respected generally. The following hierarchy pattern for the polarization fractions seems compatible with the PQCD predictions and the current data, although measurements for the B (s) → K * + K * − decays are not yet available.

C.
Triple product asymmetries The TPAs involve the perpendicular polarization amplitudes A ⊥ , which are power-suppressed relative to the longitudinal ones in the naive expectation. As stated in the previous subsection, this suppression is not numerically realized in penguin-dominated modes, so the corresponding TPAs may be notable. The predicted TPAs for the B 0 (s) → (Kπ)(Kπ) decays collected in Table IV imply that the values for most of the TPAs can reach 10 −2 in magnitude. The smallness of A 2 T is attributed to the suppression from the strong phase difference between the perpendicular and parallel polarization amplitudes, which was found to diminish in the PQCD framework [31,32]. Hence, observations of A 2 T with large values would signal new physics beyond the SM. Because the decay B 0 → (K 0 π + )(K 0 π − ) receives a contribution solely from the weak annihilation, its tiny transverse polarization component (see Table I) leads to the negligible TPAs as exhibited in Table IV. The pure penguin decays B 0 (s) → (K + π − )(K − π + ) depend on a single weak phase effectively, so A i T (true) are predicted to be zero under CP conservation. Namely,Ā i T = −A i T holds, and A i T are fake asymmetries actually. The searches for the true TPAs in the B 0 s → (K + π − )(K − π + ) decay [26] in the untagged data sample indeed show no manifest deviation from zero. The diminishing A 2 T (true) for the B 0 s → (K 0 π + )(K 0 π − ) mode can be also understood through the similarity between the perpendicular and parallel polarization amplitudes. The predicted nonvanishing fake TPAs are due to the effect of strong interactions [12], which can be tested if flavor tagged measurements are available in the future.

D. Direct CP asymmetries and S-wave-induced direct CP asymmetries
As stated before, no direct CP asymmetries are generated in the pure penguin decays B 0 (s) → (K + π − )(K − π + ) in the PQCD approach, for their decay amplitudes are governed by a single product of the CKM matrix elements. The up and charm quark penguins may contribute a weak phase difference and a direct CP asymmetry, but this effect is negligible and not taken into account here [118]. The other three modes receive the addition tree contributions induced by the operators O 1,2 , so direct CP asymmetries arise from the interference between the tree and penguin amplitudes. The numerical results of the direct CP asymmetries A dir h are listed in Tables V, where we have combined the theoretical uncertainties by adding them in quadrature as in Table IV.
The factorizable emission amplitude F LL e dominates the B 0 s → (K 0 π + )(K 0 π − ) mode, which exists in both tree and penguin contributions as shown in Eq. (45). Although the tree amplitude is CKM suppressed relative to the penguin one, the large Wilson coefficient C 2 + C 1 /3 compensates this suppression partly, such that the interference between the tree and penguin contributions may induce sizable direct CP asymmetries. It is clear from Table V that the involved components have larger asymmetries except A SS and A V S , whose values are predicted to be about one  percent. For the B 0 and B + meson decays, the tree operators O 1,2 contribute only via the annihilation amplitudes, which are power-suppressed compared with the emission ones. The CKM suppression of the tree contribution is not effective here, because |V * ub V ud | and |V * tb V td | are of similar order of magnitude, ie., O(λ 3 ). The B 0 meson decay occurs only through the annihilation diagrams, so the tree contribution, despite of being color-suppressed, is comparable with the penguin one. On the contrary, the tree amplitude in the B + mode is color-favored, and the factorizable emission diagram also enhances the penguin contribution. Consequently, large direct CP asymmetries are expected in some components of these two modes as seen in Table V. Next we compare the direct CP asymmetries extracted for the two-body decays with the previous studies based on the two-body formalism. The numerical outcomes are summarized in Table VI, which shows many differences among the predictions in various approaches. Our value for the B 0 s → K * + K * − channel is greater than the previous PQCD ones, but closer to those from the QCDF, SCET and FAT approaches. It is found that all the predicted direct CP asymmetries have the same sign, though some of the errors are still large. The two distinct PQCD results in Refs. [31] and [118] for the direct CP asymmetry A dir CP (B 0 → K * + K * − ) are attributed to the different models of the vector meson DAs. Our prediction is located in between, and has a minus sign the same as of the latter. The QCDF analysis [120] shows that the strong phases for the helicity amplitudes with h = 0, , ⊥ are either 0 or π, so the direct CP asymmetry vanishes. A similar observation was made in the same framework with the endpoint parameters being fixed by the best fit [126]. The time-like penguin annihilation contribution was neglected in the FAT approach [37] due to the lack of enough data for its determination. Because the W -exchange diagrams were neglected, the FAT [37] approach yielded A dir CP (B 0 → K * + K * − ) = 0. The SCET prediction for this quantity is absent, for the annihilation diagrams, counted as a next-to-leading-power contribution [36], were dropped. We predict a negative A dir CP (B + → K * +K * 0 ), which matches those from PQCD [118] and FAT [37] in sign, but has a lower magnitude. However, the PQCD [31], QCDF [120], and SCET [36] results have a positive sign, with the magnitudes ranging from 9.5% to 23.0%. These different predictions can be discriminated by more precise measurements in the future. It is worth stressing that the rather distinct A dir CP in our and previous PQCD studies may suggest a stronger sensitivity of direct CP asymmetries to the widths of resonant states.  [120] 9.5 ± 10.6 −24.8 ± 2.6 -15 [118] The predictions for the asymmetries from the interference with the amplitude A S + defined in Eq. (59) and for the S-wave-induced direct CP asymmetries are presented in Table VII. Because the fractions f S + in Table II are less than 10%, most values of A i D andĀ i D are only a few percent. Their constructive interference may result in slightly larger CP asymmetries, such that A 1 S in the B 0 s → (K 0 π + )(K 0 π − ) decay and A 4 S in the B 0 → (K 0 π + )(K 0 π − ) decay are over 15% in magnitude. It is encouraged to search for the S-wave-induced direct CP asymmetries in these two modes. The pure penguin decays B 0 (s) → (K + π − )(K − π + ) have vanishing A i S for a reason the same as for their vanishing A i T (true) in Table IV. The difference between A 1 S and A 3 S is ascribed to the distinct combinations of Re[A S + A * 0 ] and Re[A S + A * SS ] in Eq. (59). Table I shows that the parallel component of the pure annihilation mode B 0 → (K 0 π + )(K 0 π − ) is only of order 10 −10 , explaining why A 2 D ,Ā 2 D and A 2 S are approximately zero. Since the S-wave-induced direct CP asymmetries in the considered decays still acquire less theoretical and experimental attention, we will wait for the confrontation with future data. VII: PQCD predictions for the asymmetries (%) from the interference with the amplitude A S + and for the S-waveinduced direct CP asymmetries (%) of the considered decays in the Kπ invariant mass window of 150 MeV around the K * (892) resonance.

VI. CONCLUSION
In this paper the four-body decays B (s) → (Kπ) S/P (Kπ) S/P have been investigated under the quasi-two-body approximation in the perturbative QCD framework. We have focused on the physical observables derived in the Kπ invariant mass window of 150 MeV around the K * (892) resonance. The strong dynamics associated with the hadronization of the Kπ pairs was parametrized into the nonperturbative Kπ two-meson DAs, which include both resonant and nonresonant contributions. The evaluation of the hard kernels then reduces to the one for two-body decays. The above simplified formalism, appropriate for the leading-power regions of phase space, has been applied to three-body B meson decays successfully, and extended to four-body decays here for the first time. Because the Kπ pair can be produced in the S-wave configuration, the interference between the S-and P -wave contributions stimulate rich phenomenology in these four-body B meson decays. To describe the interference accurately, we have introduced the meson and parton fractional momenta with the final-state meson mass dependence. The relations between the final-state meson momentum fractions in the B meson rest frame and the helicity angles in the Kπ pair rest frame were given, which facilitate the discussion of various asymmetries in angular distributions.
We have computed six helicity amplitudes allowed by angular momentum conservation, and the corresponding CP -averaged branching ratios. The branching ratios of the B s modes are of order 10 −6 , and those of the B modes are at least lower by an order of magnitude because of the smaller CKM matrix elements. The single and double S-wave contributions were found to be substantial in the chosen invariant mass region and consistent with the LHCb data. We have extracted the two-body B → K * K * branching ratios from the results for the corresponding four-body decays by applying the narrow-width approximation. The obtained two-body branching ratios, except B(B 0 → K * 0K * 0 ), are slightly above the previous derivations in the two-body formalism within theoretical uncertainties, and agree with the measured values. The predicted hierarchy pattern for the longitudinal polarization fractions in the B (s) meson decays into the P -wave Kπ pairs is compatible with the data roughly. However, the different polarization data between the B s → K * 0K * 0 and B 0 → K * 0K * 0 modes, which are related to each other by the U -spin symmetry, remain puzzling, and call for more in-depth explorations.
We have also calculated the TPAs, direct CP asymmetries, and S-wave induced direct CP asymmetries in the B (s) → (Kπ) S/P (Kπ) S/P decays. It was observed that the true TPAs in most of the considered channels are tiny, of order 10 −2 or even lower. The magnitudes of the fake TPAs could be larger, but they do not reflect CP violation in B meson decays. The direct CP asymmetries in some helicity configurations, such as those of the B 0 → (K 0 π + )(K 0 π − ) and B + → (K 0 π + )(K + π − ) modes, were found to be sizable. The direct CP asymmetries in the two-body decays were also extracted from the results for the corresponding four-body ones, some of which differ from the previous PQCD analyses in the two-body framework, indicating a strong sensitivity of direct CP asymmetries to the widths of resonant states. At last, we observed significant S-wave-induced direct CP asymmetries in the pure annihilation decay B 0 → (K 0 π + )(K 0 π − ).
In our numerical study the representative theoretical uncertainties from the hadronic parameters in B meson DAs, the Gegenbauer moments of the S-and P -wave Kπ two-meson DAs, and the hard scales were taken into account, among which the hard scales are identified as the major source. It suggests that higher-order contributions to fourbody B meson decays need to be included to reduce the sensitivity to the variation of these scales. At the same time, the Gegenbauer moments should be further constrained to improve the precision of theoretical predictions, and to sharpen their confrontation with future data.