Study of $B_{(s)}^0 \to \phi\phi \to (K^+K^-)(K^+K^-)$ decays in the perturbative QCD approach

In this work, we make a detailed analysis on the penguin-dominant processes $B_{(s)}^0 \to \phi\phi \to (K^+K^-)(K^+K^-)$ in the perturbative QCD (PQCD) approach. In addition to the dominant $P$-wave resonance, the scalar background $f_0(980) \to K^+K^-$ is also accounted for. We improve the Gegenbauer moments in $KK$ two-meson distribution amplitudes by fitting the PQCD factorization formulas to measured branching ratios of three-body and four-body $B$ decays. We extract the branching ratios of two-body $B_{(s)}^0 \to \phi\phi$ decays from the corresponding four-body decay modes and calculate the relevant polarization fractions together with two relative phases $\phi_{\parallel,\perp}$, which are consistent with the previous theoretical predictions. The PQCD predictions for the"true"triple product asymmetries (TPAs) are zero which are expected in the standard model due to the vanishing weak phase difference, and support the current data reported by the CDF and LHCb Collaborations. A large"fake"TPA $\mathcal{A}_\text{T-fake}^1=30.4\%$ of the decay $B^0_s \to \phi\phi \to (K^+K^-)(K^+K^-)$ is predicted for the first time, which indicates the presence of the significant final-state interactions. The TPAs of the rare decay channel $B^0 \to \phi\phi\to (K^+K^-)(K^+K^-)$ are also predicted and can be tested in the near future.


I. INTRODUCTION
In the standard model (SM), studies of the polarization amplitudes and triple product asymmetries in the flavour-changing neutral current decays provide powerful tests for the presence of physics beyond the SM [1][2][3][4][5][6][7][8][9], especially for the decay B 0 s → φφ via a b → sss penguin process, where the φ(1020) is implied throughout the remainder paper.The B 0 s → φφ decay is a pseudoscalar to vector-vector transition, where φ is reconstructed in the K + K − final states.According to the angular momentum conservation, there are three possible spin configurations corresponding to the polarizations of the final-state vector mesons: longitudinal polarization (A 0 ), and transverse polarization with spins parallel (A ) or perpendicular (A ⊥ ) to each other.The first two states A 0 and A are CP even, while the last one A ⊥ is CP odd.Polarization amplitudes can be measured by analyzing angular distributions of final-state particles.In the factorization assumption, the longitudinal polarization should dominate based on the quark helicity analysis [10,11].In sharp contrast to these expectations, large transverse polarization of order 50% is observed in B → K * φ, B → K * ρ and B s → φφ decays [12][13][14][15][16], which poses an interesting challenge for the theory.
Interference between the CP -even (A 0 , A ) and CP -odd (A ⊥ ) amplitudes can generate asymmetries in angular distributions, the triple product asymmetries, which may signal unexpected CP violation due to physics beyond the SM.In recent years, TPAs have already been measured by Belle, BABAR, CDF and LHCb [12,13,[17][18][19][20][21][22][23][24].These triple products are odd under the time reversal transformation (T ), and also constitute potential signals of CP violation due to the CP T theorem.As we know, a non-vanishing direct CP violation needs the interference of at least two amplitudes with a weak phase difference ∆φ and a strong phase difference ∆δ.The direct CP violation is proportional to sin ∆φ sin ∆δ, while TPAs go as sin ∆φ cos ∆δ.
The key point is that the direct CP violation can only be produced if there is a nonzero strong phase difference.It has been argued that all strong phases in B decays should be rather small due to the fact that the b-quark is heavy [3].Hence, if the strong phases are quite small, the magnitude of the direct CP violation is close to zero, but the TPA is maximal.It implies that direct CP violation and TPAs complement each other.Since no tree level operators can contribute to four-body decays B 0 (s) → φφ → (K + K − )(K + K − ), there is no direct CP violation in such decay modes.However, T -odd triple products (also φ φ ϕ FIG.1: Graphical definitions of the helicity angles θ1, θ2 and ϕ for the B 0 s → φφ decay, with each quasi-two-body intermediate resonance decaying to two pseudoscalars (φ → K + K − ).θ1,2 is denoted as the angle between the direction of motion of K − in the φ rest frame and φ in the B 0 s rest frame, and ϕ is the angle between the two planes defined by K + K − in the B 0 s rest frame.
called "fake" TPAs), which are proportional to cos ∆φ sin ∆δ, can provide useful complementary information.Thus, it may be more promising to search for TPAs than direct CP asymmetries in b → s penguin decays.B 0 (s) → φφ decays are usually treated as two-body final states on the theoretical sides and have been studied in the two-body framework using QCD factorization (QCDF) [25][26][27], the perturbative QCD approaches [28][29][30], the soft-collinear effective theory (SCET) [31] and the factorization-assisted topological amplitude approach (FAT) [32].While they are at least four-body decays on the experimental side shown in Fig. 1, since the vector meson φ decays via the strong interaction with a nonzero width.The four-body B meson decays are indeed more challenging than two-body decays, but provide a number of theoretical and phenomenological advantages.On the one hand, the four-body decay amplitudes depend on five kinematic variables: three helicity angles and two invariant masses of final meson pairs, while the kinematics of two-body decays is fixed.On the other hand, the four-body decays not only receive the resonant and nonresonant contributions, but also involve the possible significant final-state interactions (FSIs) [33][34][35].
Four-body decays are still mostly unexplored from the theoretical point of view since the factorization formalism that describes a multi-body decay in full phase space is not yet available at present.Recent studies on three-body hadronic decays of B mesons based on the symmetry principles [36][37][38][39][40][41], the QCDF [42][43][44][45][46][47][48][49][50][51][52] and the PQCD approaches [53][54][55][56][57][58][59][60][61][62][63][64][65][66][67][68][69] look promising.It has been proposed that the factorization theorem of three-body B decays is approximately valid when two particles move collinearly and the bachelor particle recoils back [70,71].More details can also be found in Refs.[72,73].This situation exists particularly in the low ππ or Kπ invariant mass region ( 2 GeV) of the Dalitz plot where most resonant structures are seen.The Dalitz plot is typically dominated by resonant quasi-two-body contributions along the edge.This proposal provides a theoretical framework for studies of resonant contributions based on the quasi-two-body-decay mechanism.Recently, the localized CP violation and branching fraction of the four-body decay B0 → K − π + π + π − have been calculated by employing a quasi-two-body QCDF approach in Refs.[74,75].In our previous works [76][77][78], the PQCD factorization formalism based on the quasi-two-bodydecay mechanism for four-body B meson decays have been well established.Within the framework of PQCD approach, the branching ratios and direct CP asymmetries of four-body decays B 0 (s) → ππππ have also been studied [79].As a first step, we can only restrict ourselves to the specific kinematical configurations in which each two particles move collinearly and two pairs of final state particles recoil back in the rest frame of the B meson, see Fig. 1.Naturally the dynamics associated with the pair of final state mesons can be factorized into a two-meson distribution amplitude (DA) Φ h1h2 [80][81][82][83][84][85][86].Thereby, the typical PQCD factorization formula for the considered four-body decay amplitude can be described as the form of, where Φ B is the universal wave function of the B meson and absorbs the non-perturbative dynamics in the process.The Φ KK is the two-hadron DA, which involves the resonant and nonresonant interactions between the two moving collinearly mesons.The hard kernel H describes the dynamics of the strong and electroweak interactions in four-body hadronic decays in a similar way as the one for the corresponding two-body decays.
In this work, we study the four-body decays B 0 (s) → (K + K − )(K + K − ) in the PQCD approach based on k T factorization with the relevant Feynman diagrams illustrated in Fig. 2. The invariant mass of the K + K − pair is restricted to be within ±30MeV of the known mass of the φ meson for comparison with the LHCb data [13].The effect of identical particles has been considered in our work.In the considered (K + K − ) invariant-mass range, the vector resonance φ is expected to contribute, together with the scalar resonance f 0 (980).The S and P -wave contributions are parametrized into the corresponding timelike form factors involved in the two-meson DAs.We perform a global fit of the Gegenbauer moments in two-kaon DAs associated with both longitudinal and transverse polarizations to measured branching ratios in three-body and four-body charmless hadronic B meson decays, which will be expressed in detail in the following section.With the improved two-kaon DAs, we calculate the branching ratios and polarization fractions of each partial waves.In addition, triple-product asymmetries corresponding to the interference of the CP -odd amplitudes with the other CP -even amplitudes are predicted.The rest of the paper is organized as follows.The kinematic variables for four-body hadronic B meson decays are defined in Sec.II.The considered S and P -wave two-meson DAs are also parametrized, whose normalization form factors are assumed to take the Flatté and relativistic Breit-Wigner (BW) models [87,88].We explain how to perform the global fit, present and discuss the numerical results in Sec.III, which is followed by the Conclusion.The Appendix collects the explicit PQCD factorization formulas for all the decay amplitudes.

A. Kinematics
Considering the four-body decay B(p ), as usual, we will work in the B meson rest frame.By employing the light-cone coordinates, we define the B meson momentum p B , the total momenta of the two kaon-kaon pairs, p = p 1 + p 2 , q = p 3 + p 4 , and the quark momentum k i (i = B, p, q) in each meson in the following form: with the B meson mass m B , the parton momentum fractions x i , and the parton transverse momenta k iT , i = B, 1, 2. The explicit expressions of f ± , g ± related to the invariant masses of the meson pairs via p 2 = ω 2 1 and q 2 = ω 2 2 can be written as where . For the P -wave KK 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.The individual momenta p i (i = 1 − 4) of the four final states can be derived from the relations p = p 1 + p 2 and q = p 3 + p 4 , together with the on-shell conditions p 2 i = m 2 i for the final state meson P i or Q i , with the factors and the mass ratios r i = m 2 i /m 2 B , m i being the masses of the final state mesons.Comparing Eqs. ( 5) and ( 2), one can see that the meson momentum fractions are modified by the final state meson masses, The relation between ζ 1,2 and the polar angle θ 1,2 in the dimeson rest frame in Fig. 1 can be obtained easily, with the upper and lower limits of ζ 1,2

B. Distribution amplitudes
Without the endpoint singularities in the evaluations, the distribution amplitudes are one of the most significant nonperturbative inputs in the PQCD approach.In this section, we will briefly introduce the B meson DAs, the S-, P -wave two-kaon DAs, as well as the time-like form factors used in our calculations.In what follows the subscripts S, P are always related to the corresponding partial waves.
The light-cone hadronic matrix element for a B meson is parametrized as [89-94] where q represents a d or s quark.The two wave functions φ B and φB in the above decomposition, related to φ + B and φ − B defined in the literature [95] , obey the normalization conditions It has been shown that the contribution from φB is of next-to-leading power and numerically suppressed [90,91,96], compared to the leading-power contribution from φ B .Taking the PQCD evaluation of the B → π transition form factor F B→π 0 in Ref. [96] as an example, we find that the φB contribution to F B→π 0 is about 20% of the φ B one.The higher-twist B meson DAs have been systematically investigated in the heavy quark effective theory [97], which are decomposed according to definite twist and conformal spin assignments up to twist 6.In principle, all the next-to-leading-power sources should be included for a consistent and complete analysis, which, however, goes beyond the scope of the present formalism.Therefore, we focus only on the leading-power component with the impact parameter b being conjugate to the parton transverse momentum k BT .The B meson DA φ B (x, b) is chosen as the model form widely adopted in the PQCD approach [89][90][91][92][93]98], where the constant N B is related to the B meson decay constant f B through the normalization condition The S-wave two-kaon DA can be written in the following form [102], in which the asymptotic forms of the individual twist-2 and twist-3 components φ 0 S , φ s,t S are parametrized as [80-83] with the time-like scalar form factor F S (ω 2 ).The Gegenbauer moment a S in Eq. ( 15) is adopted the same value as that determined in Ref. [103]: a S = 0.80 ± 0.16.
The corresponding P -wave two-kaon DAs related to both longitudinal and transverse polarizations are decomposed, up to the twist 3, into [57]: The various twists φ i P in the above equations can be expanded in terms of the Gegenbauer polynomials: with the Gegenbauer coefficients a 0,T 2φ and the two P -wave form factors F P (ω 2 ) and F ⊥ P (ω 2 ).The moment a 0 2φ in the longitudinal twist-2 component φ 0 P has already been determined in a recent global analysis from the three-body B decays in the PQCD approach [67] .We will update the fitting result in the following section by taking the additional four-body decay The elastic rescattering effects in a final-state meson pair can be absorbed into the time-like form factors F ,⊥ P (ω 2 ) in the two-meson DAs according to the Watson theorem [104].For the narrow resonance φ, we usually employ the relativistic BW line shape for the form factors F P (ω 2 ) [105].The explicit formula is expressed as [88] where m φ = 1.0195GeV is the φ meson mass.The mass-dependent width Γ φ (ω) is defined as with the natural width of the φ meson Γ φ = 4.25 MeV [106].The k(ω) is the momentum vector of the resonance decay product measured in the resonance rest frame, while k(m φ ) is the value of k(ω) when ω = m φ .The explicit expression of kinematic variables k(ω) is defined in the h 1 h 2 center-of-mass frame with the Källén function λ(a, b, c) = a 2 + b 2 + c 2 − 2(ab + ac + bc) and m h1,h2 being the final state mass.The orbital angular momentum L R in the two-meson system is set to L R = 1 for a P -wave state.Due to the limited studies on the form factor F ⊥ P (ω 2 ), we use the two decay constants f (T ) φ of the intermediate particle to determine the ratio For scalar resonance f 0 (980), we adopt the Flatté parametrization where the resulting line shape is above and below the threshold of the intermediate particle [87].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.Thus, we employ the modified Flatté model suggested by D.V. Bugg [107] following the LHCb collaboration [108,109], The coupling constants g ππ = 0.167 GeV and g KK = 3.47g ππ [108,109] describe the f 0 (980) decay into the final states π + π − and K + K − , respectively.The exponential factor F KK = e −αq 2 K is introduced above the K K threshold to reduce the ρ KK factor as invariant mass increases, where q k is the momentum of the kaon in the K K rest frame and α = 2.0 ± 0.25 GeV −2 [107,108].The phase space factors ρ ππ and ρ KK read as [87,108,110]

C. Helicity amplitudes
The differential branching fraction for the meson rest frame is expressed 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 the K + K − pair in the B (s) meson rest frame.The four-body phase space has been derived in the analysis of the K → ππlν decay [111], the semileptonic B → D(D * )πlν decays [112], semileptonic baryonic decays [113,114], and four-body baryonic decays [115].One can confirm that Eq. ( 31) is equivalent to those in Refs.[113,115] by appropriate variable changes.Replacing the helicity angle θ by the meson momentum fraction ζ via Eq.( 8), the Eq. ( 31) is turned into The B 0 s → φφ → (K + K − )(K + K − ) decay comprises a mixture of CP eigenstates and can be disentangled by means of an angular analysis in the helicity basis.In this basis, the decay is described by three angles θ 1 , θ 2 and ϕ, depicted in Fig. 1, where the θ 1,2 is the angle between the K − direction in the φ rest frame and the φ direction in the B 0 s rest frame, and ϕ is the angle between the two φ meson decay planes.
Due to the proximity of the φ resonance to the scalar f 0 (980) resonance, there are irreducible scalar resonant contributions to four-body B 0 (s) → (K + K − )(K + K − ) decays.Thereby, a K + K − pair can be produced in the S or P -wave configuration in the selected invariant mass regions.One decomposes the decay amplitudes into six helicity components: h = V V (3), V S (2), and SS, each with a corresponding amplitude A h , where V denotes a vector meson and S denotes a scalar meson.The first three, commonly referred to as the P -wave amplitudes, are associated with the final states, where both K + K − pairs come from intermediate vector mesons.In the transversity basis, a P -wave decay amplitude can be decomposed into three components: A 0 , for which the polarizations of the final-state vector mesons are longitudinal 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 + K − pair can arise from R 1 or R 2 labelled in Fig. 2(a), the corresponding single S-wave amplitude is denoted A V S .The double S-wave amplitude A SS is associated with the final state, where both two-meson pairs are generated in the S wave.A randomised choice is made for which φ meson is used to determine θ 1 and which is used to determine θ 2 .Thus, the total decay amplitude A is a coherent sum of the P -, S-, and double S-wave components.Specifically, these helicity amplitudes for the B 0 (s) → (K By including the ζ 1,2 dependencies instead of θ 1,2 and azimuth-angle dependencies relying on Eq. ( 8), the total decay amplitude in Eq. ( 33) can be written as On basis of Eq. ( 33), we can obtain the branching ratio form, where the invariant masses ω 1,2 are integrated over the chosen K + K − mass window.The coefficients C h are the results of the integrations over ζ 1 , ζ 2 , ϕ in terms of Eq. ( 36) and listed as follows, The CP -averaged branching ratio and the direct CP asymmetry in each component are defined as below, respectively, where Bh 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.
For the V V decays, the polarization fractions f λ with λ = 0, , and ⊥ and two relative phases φ , φ ⊥ are described as with the normalisation relation

D. Triple product asymmetries
Consider a four-body decay B → R 1 (→ , in which one measures the four particles' momenta in the B rest frame.We define nRi (i = 1, 2) is a unit vector perpendicular to the R i decay plane and ẑR1 is a unit vector in the direction of R 1 in the B rest frame.Thus we have implying a T -odd scalar triple product One can define a TPA as an asymmetry between the number of decays involving positive and negative values of sin ϕ or sin 2ϕ, It has been found that TPAs originate from the interference of the CP -odd amplitudes A ⊥ with the other CP -even amplitudes A 0 , A .According to Eq. ( 8), the TPAs associated with A ⊥ for the considered four-body decays are derived from the partially integrated differential decay rates as [4,23] with the denominator The above TPAs contain the integrands Im(A ⊥ A * 0, ) = |A ⊥ ||A * 0, | sin(∆φ + ∆δ), where ∆φ and ∆δ denote the weak and strong phase differences between the amplitudes A ⊥ and A 0, , respectively.As already noted, Im(A ⊥ A * 0, ) can be nonzero even if the weak phases vanish.Thus, it is not quite accurate to identify a nonzero TPA as a signal of CP violation.To obtain a true CP violation signal, one has to compare the TPAs in the B and B meson decays.The helicity amplitude for the CP -conjugated process can be inferred from Eq. ( 35) through A 0 → Ā0 , A → Ā and A ⊥ → − Ā⊥ , in which the Āλ are obtained from the A λ by changing the sign of the weak phases.Thus, the TPAs Āi T for the charge-conjugate process are defined similarly, but with a multiplicative minus sign.[28,67].Other parameters are from PDG 2020 [106].One therefore constructs the "true" and "fake" asymmetries by combining A i T and Āi T [4] with Γ being the decay rate of the CP -conjugate process, T = (2ζ 1 − 1)(2ζ 2 − 1) sin ϕ and T = (2ζ 1 − 1)(2ζ 2 − 1) sin φ and the denominator is for the CP -conjugate decay.
It is shown that the terms Im[A ⊥ A * 0, − Ā⊥ Ā * 0, ] in Eqs. ( 48) and ( 49) are proportional to sin ∆φ cos ∆δ, which are nonzero only in the presence of the weak phase difference.Then TPAs provide 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 is maximal when the strong phase difference vanishes.For the special case of the involved neutral intermediate states B 0 (s) → φφ modes, in which each helicity amplitude involves the same single weak phase in the SM.This results in A i T-true = 0 due to the vanishing weak phase difference.The "true" TPAs for the neutral modes are thus expected to be zero in the SM.If such asymmetries are observed experimentally, it is probably a signal of new physics.While for the term Im[A ⊥ A * 0, + Ā⊥ Ā * 0, ] ∝ cos ∆φ sin ∆δ, the A i T-fake can be nonzero when the weak phase difference vanishes.Such a quantity is referred as a fake asymmetry (CP conserving), which reflects the effect of strong phases [4,5], instead of CP violation.

III. NUMERICAL ANALYSIS
In this section, we calculate the branching rations (B), the polarization fractions f λ and relative phases φ ,⊥ (rad), together with TPAs, respectively.The related input parameters for the numerical calculations are collected in Table I.The decay constants are used the values from Refs.[28,67], while the meson masses, Wolfenstein parameters, and the lifetimes are taken from the PDG review [106].We neglect uncertainties on the constants since they are negligible with respect to other sources of uncertainties.

A. Global fit
According to Eqs. ( 20)-( 25), the total amplitudes A related to both longitudinal (L) and transverse (N, T ) components for the four-body decays B 0 (s) → φφ → (K + K − )(K + K − ) can be expanded in terms of the Gegenbauer moments from the two-meson DAs.As a result, we can decompose the squared amplitudes into the linear combinations of the Gegenbauer moments a 0,T 2φ and their products While for three-body decays B (s) → (π, K)φ → (π, K)KK, analogously, the squared amplitudes |A| 2 can be paramatrised as follows We then compute the coefficients M , which involve only the Gegenbauer polynomials, to establish the database for our global fit.
Similar to the proposal in Refs.[67,99], we adopt the standard nonlinear least-χ 2 (lsq) method [116], in which the χ 2 function is defined for n pieces of experimental data v i ± δv i with the errors δv i and the fitted corresponding theoretical values v th i as In general, we should include maximal amount of data in the fit in order to minimize statistical uncertainties.However, those measurements with significance lower than 3σ do not impose stringent constraints, and need not be taken into account in principle.Therefore, the Gegenbauer moments a

0(T ) 2φ
for the twist-2 KK DAs φ 0(T ) P can be obtained by fitting the formulas in Eqs. ( 53)-( 55) with the Gegenbauer-moment-independent database to the five pieces of B 0(+) → K 0(+) φ → K 0(+) KK and B 0 s → φφ → (K + K − )(K + K − ) data, including three branching ratios and two polarization fractions as summarized in Table II, whose errors mainly arise from experimental uncertainties.For comparison, the updated fitting results are also listed in Table II and match well with the data within errors.Note that our a 0 2φ , determined with χ 2 /d.o.f.= 1.3, is distinct from the value a 0 2φ = −0.31± 0.19 in Ref. [67], which can be understood from the following clarification.The additional new four-body decay B 0 s → φφ → (KK)(KK) included in the present work is dominated by B 0 s → (φ →)KK transition form factors, the B of which could be more sensitive to the Gegenbauer moment a 0 2φ .Hence, the measured B 0 s → φφ → (KK)(KK) branching ratio can give an effective constraint on the global fit of the KK two-meson DAs, and the corresponding fitting result of the a 0 2φ could be changed a lot: from −0.31 ± 0.19 to 0.40 ± 0.06.
One can also observe that the a T 2φ fitted in this work is slightly larger than unity as shown in Eq. ( 57), which is not favored in view of the convergence of the Gegenbauer expansion.We then added one more Gegenbauer moment a T 4φ in the twist-2 transverse component φ T P .Naturally, the Eq. ( 54) should be replaced with the following form, and a fit with χ 2 /d.o.f.= 1.3 is attained, The outcome of both a T 2φ and a T 4φ in Eq. ( 59) are all smaller than unit, implying that the contributions from the higher-order Gegenbauer moment is significant.In principle, we should introduce the same number of Gegenbauer moments for the two twist-2 DAs φ 0 P and φ T P .However, it is not practical to include many parameters in the fit because of the limited amount of experimental data at present.Consequently, we will adopt the a 0,T 2φ as presented in Eq. ( 57) in this work.Anyway, the above results show that the a T 2φ can be reduced efficiently by including the higher moment a T 4φ into the fit when more experimental data with improved precision are available in the future.

B. S-wave contributions
TABLE III: PQCD predictions for the branching ratios of various components and their sum in the B 0 (s) → (K + K − )(K + K − ) decays.The theoretical uncertainties are attributed to the variations of the shape parameter ωB (s) in the B (s) meson DA, of the Gegenbauer moments in various twist DAs of KK pair, and of the hard scale t and the QCD scale ΛQCD.
(1.73   The PQCD predictions for the branching ratios of various components and their sum in the B 0 (s) → (K + K − )(K + K − ) decays are summarized in Table III, in which the theoretical uncertainties are estimated from three different sources.The first error is due to the shape parameters ω B in the B (s) meson DAs with 10% variation.The second one comes from the Gegenbauer moments in various twist DAs of KK pair with different intermediate resonances.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 three uncertainties are comparable, and their combined impacts could exceed 50%, implying that the nonperturbative parameters in the DAs of the initial and final states need to be constrained more precisely, and the higher-order correction to four-body B meson decays is critical.It should be stressed that these considered modes are induced only by penguin operators in the PQCD approach at leading order as can be seen easily from the Appendix, their direct CP violations are naturally zero without the interference between the tree and penguin amplitudes.
In contrast to the vector mesons, the identification of scalar mesons is a long-standing puzzle, and the underlying structure of scalar mesons is not theoretically well established (for a review, see Ref. [106]).Based on the assumption that f 0 (980) is a pure ss state, different kinds of theoretical approaches have been applied to study the B (s) meson decays involving f 0 (980) in the final states, for instance: (a) the charmonium decay B s → J/ψf 0 (980) was analysed in the light-cone QCD sum rule and factorization assumption [117], as well as the generalized factorization and SU (3) flavor symmetry [118]; (b) in Ref. [119], the Bs → f 0 (980) transition form factor was calculated from the light-cone sum rules with B-meson DAs.As a first approximation, the scalar meson f 0 (980) is taken into account in the ss density operator in the present work.The S-wave time-like form factor F S (ω 2 ) adopted to parameterize the S-wave two-kaon DAs have been determined in Refs.[66,78,103].
By using the numerical results given in Table III, the S-wave fractions defined as are also calculated in our work, and the corresponding results are presented in Table IV.The total S-wave fraction of the two considered modes is estimated to be less than 2%, which is in good agreement with the experimental analysis that the contributions from the S-wave components are negligible [13].
The double S-wave decays B 0 (s) → f 0 (980)f 0 (980) have already been systematically studied in the two-body framework within the PQCD approach [127,128].Taking the B 0 s → f 0 (980)f 0 (980) decay as an example, we can roughly estimate III on basis of Eq. ( 61).It is worthwhile to note that the estimation B = (1.43 +1.29 −0.91 ) × 10 −7 is much smaller than the previous two PQCD results B(B 0 s → f 0 (980)f 0 (980)) = (2.66 +1.08 −0.85 ) × 10 −4 [127] and B(B 0 s → f 0 (980)f 0 (980)) = (5.31+1. 74 −1.39 ) × 10 −4 [128].Strictly speaking, the narrow width approximation is actually not fully justified since such approximation has its scope of application.As pointed out in Refs.[129,130], the narrow width approximation should be corrected by including finite-width effects for the broad scalar intermediate state.The two-body result extracted from the four-body branching ratio may suffer from a large uncertainty due to the finite-width effects of the scalar resonance.In addition, as stated above, the S-wave contributions show a strong dependence on the range of the KK invariant mass.We hope that the future LHCb and Belle II experiments can perform a direct measurement on four-body decays B 0 s → φf 0 (980 C. Branching ratios and polarization fractions of two-body B 0 (s) → φφ decays With the narrow width approximation in Eq. ( 61), the branching ratios of two-body B 0 (s) → φφ decays are extracted in Table V.The polarization fractions of the two-body B 0 (s) → φφ decays together with two relative phases calculated in this work TABLE V: Branching ratios, polarization fractions and relative phases for the two-body B 0 (s) → φφ decays.For comparison, we also list the results from PQCD [28,29], QCDF [26], SCET [31], and FAT [32].The world averages of experimental data are taken from PDG 2020 [106].The sources of the theoretical errors are the same as in Table III   3.50 ± 0.17 are also listed in Table V.For a comparison, we display the updated predictions in the QCDF [26], the previous predictions in the PQCD approach [28,29], SCET [31] and FAT [32].Experimental results for branching ratios, polarization fractions and relative phases are taken from PDG 2020 [106].
It is obvious that most of the theoretical predictions of B 0 s → φφ decay are consistent well with experiments within errors.In Ref. [29], the authors kept the additional power corrections related to the ratio r 2 V = m 2 V /m 2 B (m V and m B denote the masses of the vector and B mesons, respectively).By including the r 2 V term, their result B(B 0 s → φφ) = (16.7 +4.9 −3.8 ) × 10 −6 is about twice smaller than B(B 0 s → φφ) = (35.3+18.7 −12.3 ) × 10 −6 [28], and more close to the experimental data B(B 0 s → φφ) = (18.7 ± 1.5) × 10 −6 [106].As stated in Refs.[53,54], the branching ratios in the quasi-two-body mechanism show their dependence on the invariant masses of the final-state meson pairs.In this work, the factors [29] when ω 1,2 = m φ .In addition, the partonic kinematic variables have been refined to take into account finite masses of final-state mesons, which can suppress the branching ratio of the B 0 s → φφ decay effectively.Therefore, our result agrees well with that from the updated PQCD [29].
The rare decay B 0 → φφ can occur only via penguin annihilation topology in the SM.The predicted branching ratio is very small at (10 −8 ), which makes it sensitive to any new physics contributions.The current experiment gives the upper limit: B(B 0 → φφ) < 2.7 × 10 −8 at 90% CL [13], so the more accurate experimental results are needed to test the theory.It is observed that the branching fraction of B 0 → φφ decay is much smaller than that of B 0 s → φφ decay by almost three orders.There are two main reasons: the one is that the B 0 → φφ governed by b → d transition is highly suppressed by the CKM matrix elements |V td /V ts | 2 ∼ 0.05, the other is that B 0 → φφ belongs to the pure annihilation decay.As is known, the contributions from the annihilation diagrams (Figs.2(e)-2(h)) are always power suppressed compared to the factorizable emission diagrams (Figs.2(a) and 2(b)) in the PQCD approach .What's more, there exists a big cancellation between the two factorizable annihilation diagrams Figs.2(e) and 2(f) for the contributions from the (V − A)(V − A) operators, especially when two final state mesons are identical, like B 0 → φφ.
For the charmless B decays, it is naively expected that the helicity amplitudes H i (with helicity i = 0, −, +) satisfy the hierarchy pattern which are related to the spin amplitudes (A 0 , A , A ⊥ ) in Appendix by The above hierarchy relation satisfies the expectation in the factorization assumption that the longitudinal polarization should dominate based on the quark helicity analysis [10,11].In sharp contrast to these expectations, roughly equal longitudinal and transverse components are found in measurements of B → K * φ, B → K * ρ decays [12,[14][15][16].Measurements of the low longitudinal polarization fraction in B 0 s → φφ by CDF [19] and LHCb [13,20,22] indicate a large transverse polarization.The longitudinal polarization fraction f 0 = 0.381 ± 0.007 ± 0.012 of the B 0 s → φφ decay has been reported by LHCb recently [13], where the first uncertainty is statistical and the second systematic.This shows that the scaling behavior shown in Eq. ( 64) is violated.The interest in the polarization in penguin transition, such as b → s decay B 0 s → φφ, is motivated by its potential sensitivity to physics beyond the SM.
As shown in Table V, we obtain the longitudinal polarization fraction f 0 = (38.2+7.2 −7.6 )% of the B 0 s → φφ decay with the updated Gegenbauer moments of two-meson DAs, which is consistent with the previous PQCD calculation [29] and those from QCDF [26], SCET [31] and FAT [32] within uncertainties.In the PQCD approach, the large transverse polarization fraction can be interpreted on the basis of the chirally enhanced annihilation diagrams, especially the (S − P )(S + P ) penguin annihilation, introduced by the QCD penguin operator O 6 [131], which is originally introduced in Ref. [132].A special feature of the (S − P )(S + P ) penguin annihilation operator is that the light quarks in the final states are not produced through chiral currents.So, there is no suppression to the transverse polarization caused by the helicity flip.Then the polarization fractions satisfy For the B 0 → φφ decay, the longitudinal polarization contribution is dominant, which is consistent with the recent updated PQCD calculation [29] and also verified in Ref. [133].As clarified before, the contributions from the factorizable annihilation diagrams (Figs. ) components are all power suppressed.Thus, the total transverse contributions are actually negligible, leading to f 0 ∼ 1 as shown in Table V.
The relative phases φ and φ ⊥ of the B 0 (s) → φφ decays are also studied in the present work as shown in Table V.In fact, two relative phases derived from the decay amplitudes A 0, ,⊥ in Eq. ( 40) are dependent on the invariant mass ω 1,2 .We fix ω 1,2 = m φ in our calculation for comparison with the two-body analysis.For B 0 s → φφ decay, the differences between the two previous PQCD results [28,29] are mainly attributed to the treatment of the terms in the decay amplitude proportional to the ratio V /m 2 B , which has been neglected in Ref. [28].While in our calculations, the factors η 1,2 = ω 2 1,2 /m 2 B given in Eq. ( 3) become equal to r 2 V = m 2 V /m 2 B in Ref. [29] when ω 1,2 = m φ .Therefore, our new four-body computation φ = (1.86 +0.17 −0.24 ) rad and φ ⊥ = (1.85 ± 0.18) rad agree with the updated PQCD results φ = (2.01 ± 0.23) rad and φ ⊥ = (2.00+0.24 −0.21 ) rad [29] within errors.It is obvious that our predictions of the relative phases are smaller than those of the SCET [31] and FAT [32] calculations as well as the experimental data [106].As stressed above, we use the two decay constants f (T ) φ of the intermediate particle to determine the ratio F ⊥ P (ω 2 )/F P (ω 2 ) ≈ (f T φ /f φ ).To be honest, we have omitted the phase difference between the two form factors F ⊥ P (ω 2 ) and F P (ω 2 ) due to the limited studies on the form factor F ⊥ P (ω 2 ).We have found that the gap between our predictions and the measurements of two relative phases can be resolved effectively by introducing an additional phase β in the above approximate equation, In  [29].The main reason is that we have kept track of the additional higher power corrections related to the momenta fraction x B , which has been ignored in Ref. [29].We have reexamined the two phases φ ,⊥ without the contributions from the x B : φ = 2.67 rad, φ ⊥ = 2.75 rad, which are similar to the two body analysis.It should be stressed that the contributions from the annihilation diagrams (Figs.2(e)-2(h)) are of higher power themselves for a pure annihilation decay mode without chiral enhancement, like B 0 → φφ.In that case, the terms proportional to x B in the amplitudes are not negligible and should be reserved in the calculations.Anyway, all these theoretical predictions need to be further tested in the future when more data are available.
In the involved neutral intermediate states B 0 s → φφ and B 0 → φφ modes, each helicity amplitude involves the same single weak phase in the SM.This results in A i T-true = 0 due to the vanishing weak phase difference.The "true" TPAs for these neutral modes are thus predicted to be zero in the SM as shown in Table VII.If such asymmetries are observed experimentally, it is probably signify the presence of new physics.On the experimental side, the measurements of TPAs for B 0 s → φφ → (K + K − )(K + K − ) have been reported by CDF [19] and LHCb Collaborations [20,22] and have shown no evidence of deviations from the SM.The most recent measurements of the "true" TPAs give [13] A V = −0.014± 0.011(stat) ± 0.004(syst), A U = −0.003± 0.011(stat) ± 0.004(syst), where the first uncertainty is statistical and the second systematic.No evidence for CP violation is found, which is consistent with SM predictions.The predicted "fake" TPAs for the B 0 (s) → (K + K − )(K + K − ) decays are presented in Table VII.As "fake" TPAs are due to strong phases and require no CP violation, the large fake A 1,2  T-fake simply reflects the importance of the strong final-state phases.The magnitude of A 1 T-fake for the B 0 s channel exceeds ten percent and reaches 30.4%.The sizable magnitude is mainly enhanced by the strong phase difference between the longitudinal and perpendicular polarization amplitudes, which is found in Table V.The smallness of A 2 T-fake is attributed to the suppression from the strong phase difference between the perpendicular and parallel polarization amplitudes, which can be seen in Table V and has been verified by LHCb [13].Hence, observations of A 2 T-fake with large values would signal new physics beyond the SM.As mentioned above, the hierarchy in Eq. ( 64) is numerically not respected by penguin-dominated decays.Thus, final states with large transverse amplitude fractions are favourable for the measurement of TPAs and can provide valuable complementary information on CP violation without requiring the generation of a sizable strong phase difference.Our predictions can be tested in the future.

IV. CONCLUSION
In this work, we have studied the related helicity amplitudes of four-body B 0 (s) → (K + K − )(K + K − ) decays based on the angular analysis, where K + K − invariant-mass spectrum is dominated by the vector resonance φ.The scalar resonance f 0 (980) is also contributed in the K + K − invariant-mass range.The strong dynamics of the scalar or vector resonance decays into the meson pair is parametrized into the corresponding two-meson distribution amplitude, which has been established in three-body B meson decays and further improved by performing a global fit through combining the measured branching ratios in four-body decays.
The branching ratios of four-body B 0 (s) → (K + K − )(K + K − ) decays are presented with the updated P -wave two-kaon distribution amplitudes.We have extracted the two-body B 0 (s) → φφ branching ratios from the results for the corresponding four-body decays under the narrow-width approximation and shown the polarization fractions and relative phases of the decay channels.The obtained two-body branching ratios agree well with previous theoretical studies in the two-body framework within errors.The predicted hierarchy pattern for the longitudinal polarization fractions in the B 0 (s) → φφ decays is in agreement with the data.
Since the triple product asymmetries are helpful to discover physics beyond the standard model, we perform an angular analysis and estimate the triple product asymmetries on four-body B 0 (s) → φφ → (K + K − )(K + K − ) decays.The "true" TPAs of four-body B 0 s → φφ → (K + K − )(K + K − ) decays are predicted to be zero due to the vanishing weak phase difference, which is consistent with the experiments.The prediction of "fake" TPA A 1 T-fake of B 0 s → φφ → (K + K − )(K + K − ) reaches 30% in magnitude, which reflects the importance of the strong final-state phases and can be tested in the future.We also make predictions of TPAs for B 0 → φφ → (K + K − )(K + K − ) decay and wait for the confrontation with future data.

Appendix A: Decay amplitudes
In the Appendix, we present the PQCD factorization formulas for the amplitudes of the considered four-body hadronic B meson decays: • B → f 0 (980)f 0 (980) → (K + K − )(K + K − ) decay modes where G F = 1.16639 × 10 −5 GeV −2 is the Fermi coupling constant and the V ij 's are the Cabibbo-Kobayashi-Maskawa matrix elements.The superscripts LL, LR, and SP refer to the contributions from (V − A) ⊗ (V − A), (V − A) ⊗ (V + A), and (S − P ) ⊗ (S + P ) operators, respectively.The explicit formulas for the factorizable emission (annihilation) contributions F e(a) and the nonfactorizable emission (annihilation) contributions M e(a) from Fig. 2 can be obtained easily in Ref. [76].

0 s
The shape parameter takes the values ω B = 0.40 GeV for B 0 meson and ω Bs = 0.48 GeV [89, 99-101] for B meson with 10% variation in the numerical study below.
while the moment a T 2φ associated with the transverse twist-2 component φ T P is determined in a global analysis for the first time.Since the amounts of the current experimental data are not yet enough for fixing the Gegenbauer moments in the twist-3 DAs φ s,t P and φ v,a P , they have been set to the asymptotic forms in the present work.

Mass
2(e)-2(f)) are canceled by each other because of the current conservation.Hence, the terms like F LL,0( ) aφ and F LR,0( ) aφ in the Eq.(A1) are exactly equal to zero, while the only left parts F LL,⊥ aφ and F LR,⊥ aφ for the factorizable emission diagrams are power suppressed.For the non-factorizable annihilation diagrams (Figs.2(g)-2(h)), the longitudinal parts give the leading and dominant contributions, and other terms related to parallel (M LL, aφ , M SP, aφ ) and perpendicular (M LL,⊥ aφ , M SP,⊥ aφ

TABLE I :
The decay constants are taken from Refs.

TABLE II :
(57)]imental data for branching ratios and polarization fractions[106], and the theoretical results derived from the fitted Gegenbauer moments in Eq.(57).For simplicity, only the theoretical errors from the Gegenbauer moments are presented.

TABLE IV
: PQCD predictions of the S-wave fractions in the B 0 (s) → (K + K − )(K + K − ) decays.The sources of the theoretical errors are the same as in TableIII.
but added in quadrature.

TABLE VII :
PQCD predictions for the TPAs (%) of the four-body B 0 (s) → (K + K − )(K + K − ) decays.The sources of theoretical errors are same as in TableIIIbut added in quadrature.